A Flexibly Conditional Screening Approach via a Nonparametric Quantile Partial Correlation

Abstract

Considering the influence of conditional variables is crucial to statistical modeling, ignoring this may lead to misleading results. Recently, Ma, Li and Tsai proposed the quantile partial correlation (QPC)-based screening approach that takes into account conditional variables for ultrahigh dimensional data. In this paper, we propose a nonparametric version of quantile partial correlation (NQPC), which is able to describe the influence of conditional variables on other relevant variables more flexibly and precisely. Specifically, the NQPC firstly removes the effect of conditional variables via fitting two nonparametric additive models, which differs from the conventional partial correlation that fits two parametric models, and secondly computes the QPC of the resulting residuals as NQPC. This measure is very useful in the situation where the conditional variables are highly nonlinearly correlated with both the predictors and response. Then, we employ this NQPC as the screening utility to do variable screening. A variable screening procedure based on NPQC (NQPC-SIS) is proposed. Theoretically, we prove that the NQPC-SIS enjoys the sure screening property that, with probability going to one, the selected subset can recruit all the truly important predictors under mild conditions. Finally, extensive simulations and an empirical application are carried out to demonstrate the usefulness of our proposal.

1. Introduction

Variable screening technique has been demonstrated as a computationally fast and efficient tool in solving many problems in ultrahigh dimensions. For example, in many scientific areas, such as biological genetics, finance and econometrics, we may collect the ultrahigh dimensional data sets (e.g., biomarkers, financial factors, assets and stocks), where the number $p_n$ of predictors extremely exceeds the sample size n. Theoretically, ultrahigh dimension often refers to the dimensionality $p_n$ and sample size n satisfies the relationship: $p_n=O(exp(n^a))$ for some constant $a>0$. Variable screening is able to reduce the computational cost, to avoid the instability of algorithms, and to improve the estimation accuracy. These issues exist in the variable selection approaches based on LASSO [1], SCAD [2,3] or MCP [4] for ultrahigh dimensional data. Since the seminal work of [5], which pioneeringly proposed the sure independence screening (SIS) procedure, many variable screening approaches have been consecutively documented over the last fifteen years, including the model-based methods (e.g., [6,7,8,9,10,11]) and the model-free methods [12,13,14,15,16,17,18,19,20]. These papers have showed that with probability approaching one, the set of selected predictors contain the set of all truly important predictors.

Most marginal approaches focus only on developing various effective and robust measures to characterize the marginal association between the response and individual predictor. Whereas, these methods do not take into consideration the influence of conditional variables or confounding factors on the response. A simple application of SIS is relatively rough since SIS may perform poorly when predictors are highly correlated with each other. Some predictors that are weakly relevant or irrelevant, but jointly correlated to the response, may be excluded in the final model after applying marginal screening methods. This will result in a high false positive rate (FPR). To surmount this weakness, an iterated screening algorithm or a penalization-based variable selection is usually offered as a refined follow-up step (e.g., [5,10]).

Conditional variable screening can be viewed as an important extension of the marginal screening. It accounts for conditional information when calculating the marginal screening utility. There is relatively less work in the literature. To name a few, Ref. [21] proposed a conditional SIS (CIS) procedure to improve the performance of SIS because some correlated conditional variables may increase the chance of boosting the rank of the marginally weak predictor and that of reducing the number of false negatives. The paper [22] proposed a confounder-adjusted screening method for high dimensional censoring data, in which the additional environmental confounders are regarded as conditional variables. The researchers in [23] studied the variable screening by incorporating within-subject correlation for ultrahigh dimensional longitudinal data, where they used some baseline variables as conditional variables. Ref. [24] proposed a conditional distance correlation-based screening via kernel smoothing method, while [25] further presented a screening procedure based on conditional distance correlation, which is similar to [24] in methodology, but differs in theory. Additionally, Ref. [11] developed a conditional quantile correlation-based screening approach using the B-spline smoothing technique. However, in [11,24,25], among others, the conditional variable they considered is only univariate. Further, Ref. [21] focuses on the generalized linear models, but cannot handle heavy-tailed data. For this regard, we aim to develop a screener that behaves more robustly to outliers and heavy-tailed data, and simultaneously considers more than one conditional variable. On the choice of conditional variables, one can achieve that through some prior knowledge such as published research work or the experience of experts from relevant subjects. When no prior knowledge is available, one can apply some marginal screening approaches, such as the SIS or its robust variants, to select several top-ranked predictors as conditional variables.

On the other hand, to the best of our knowledge, several works have considered multiple conditional variables based on distinct partial correlations. For instance, Ref. [26] proposed a thresholded partial correlation approach to select significant variables in linear regression models. Additionally, Ref. [17] presented a screening procedure on the basis of the quantile partial correlation in [27], and they referred to the procedure as QPC-SIS. More recently, Ref. [28] proposed a copula partial correlation-based screening approach. It is worth noting that the partial correlation used in both [17,28] removes the effect of conditional variables on the response and each predictor through fitting two parametric models with a linear structure. However, this manner may be ineffective, especially when the conditional variables have a nonlinear influence on the response nonlinear. This motivates us to work out a flexible way to control the impact of conditional variables. Meanwhile, we also take into account the issue of the robustness to outlying or heavy-tail response in this paper.

This paper contributes a robust and flexible conditional variable screening procedure via a partial correlation coefficient, which is a non-trivial extension of [17]. First of all, in order to precisely control conditional variables, we propose a nonparametric definition of QPC, which extends that of [17] and allows for more flexibility. Specifically, we first fit two nonparametric additive models to remove the effect of conditional variables on the response and an individual predictor, where we use the B-spline smoothing technique to estimate the nonparametric functions. This can be viewed as a nonparametric adjustment for controlling conditional variables. By that, we can obtain two residuals, on which a quantile correlation can be calculated to formulate a nonparametric QPC. Second, we use this quantity as the screening utility in variable screening. This procedure can be implemented rapidly. We refer to this procedure as the nonparametric quantile partial correlation-based screening, denoted as NQPC-SIS. Third, theoretically, we establish the sure screening property for NQPC-SIS under some mild conditions. Compared to [17], our approach is more flexible and our theory on the sure screening property is more difficult to derive. Moreover, our screening idea can be easily transferred to some existing screening methods that use some popular partial correlation.

The remainder of the paper is organized as follows. In Section 2, the NQPC-SIS is introduced. The technical conditions needed are listed and asymptotic properties are established in Section 3. Section 4 provides an iterative algorithm for a further refinement. Numerical studies and empirical analysis of real data set are carried out in Section 5. Concluding remarks are given in Section 6. All the proofs of the main results are relegated to the Appendix A.

2. Methodology

2.1. A Preliminary

In this section, we formally introduce the NQPC-SIS procedure. To begin with, we give some background on the quantile correlation (QC) introduced in [27]. Let X and Y be two random variables, and $EX$ be the expectation of X. The definition of QC is formulated as

$${\mathrm{qcor}}_{\tau }(Y,X)=\frac{\mathrm{E}[{\psi }_{\tau }(Y-Q_{\tau ,Y})(X-\mathrm{E}(X))]}{\sqrt{\mathrm{var}(I(Y-Q_{\tau ,Y}>0))\mathrm{var}(X)}},$$

(1)

where $Q_{\tau ,Y}$ is the $\tau$th quantile of Y, and ${\psi }_{\tau }(u)=\tau -I(u<0)$ for some quantile level $\tau \in (0,1)$, here $I(·)$ denotes an indicator function. This correlation takes on a value between $-1$ and 1 and is asymmetric with respect to Y and X compared with the conventional correlation coefficient. The QC shares the merits: the property of monotone invariance for Y as well as the robustness of Y, due to the use of the quantile rather than the mean in the definition. Thus, QC affects little in the presence of outliers in Y. Besides, as shown in [27], ${\mathrm{qcor}}_{\tau }(Y,X)$ is closely related to the quantile regression. If we denote by $(a_{0\tau }^{*},a_{1\tau }^{*})$ the minimizer of $\mathrm{E}\{{\rho }_{\tau }(Y-a_{0\tau }-a_{1\tau }X)\}$ with respect to $a_{0\tau }$ and $a_{1\tau }$, where ${\rho }_{\tau }(u)=u[\tau -I(u<0)]$. Then, it follows that ${\mathrm{qcor}}_{\tau }(Y,X)=\phi (a_{1\tau }^{*})$, where $\phi (·)$ is a continuous and increasing function, and $\phi (a_{1\tau }^{*})=0$ if and only if $a_{1\tau }^{*}=0$.

When QC is used as a marginal screening utility for variable screening, the screening results obtained may be misleading when the predictors are highly correlated. To overcome this problem, Ref. [17] proposed the screening based on quantile partial correlation (QPC) to reduce the effect from conditional predictors. For the sake of presentation, write ${\mathbf{X}}_{-j}={(X_k,k\ne j)}^T$ for $j=1,\dots ,p_n$. The QPC in [17] is defined as

$$\begin{array}{ccc}\hfill {\mathrm{qpcor}}_{\tau }(Y,X_j|{\mathbf{X}}_{-j})& =& \frac{\mathrm{cov}({\psi }_{\tau }(Y-{\mathbf{X}}_{-j}^T{\mathbf{\alpha }}_j^0),X_j-{\mathbf{X}}_{-j}^T{\mathbf{\theta }}_j^0)}{\sqrt{\mathrm{var}({\psi }_{\tau }(Y-{\mathbf{X}}_{-j}^T{\mathbf{\alpha }}_j^0))\mathrm{var}(X_j-{\mathbf{X}}_{-j}^T{\mathbf{\theta }}_j^0)}}\hfill \\ & =& \frac{\mathrm{E}\{{\psi }_{\tau }(Y-{\mathbf{X}}_{-j}^T{\mathbf{\alpha }}_j^0)(X_j-{\mathbf{X}}_{-j}^T{\mathbf{\theta }}_j^0)\}}{\sqrt{\tau (1-\tau ){\sigma }_j^2}},\hfill \end{array}$$

(2)

where ${\sigma }_j^2=\mathrm{var}(X_j-{\mathbf{X}}_{-j}^T{\mathbf{\theta }}_j^0)$, ${\mathbf{\alpha }}_j^0={\mathrm{argmin}}_{{\mathbf{\alpha }}_j}\mathrm{E}\{{\rho }_{\tau }(Y-{\mathbf{X}}_{-j}^T{\mathbf{\alpha }}_j)\}$ and ${\mathbf{\theta }}_j^0={\mathrm{argmin}}_{{\mathbf{\theta }}_j}\mathrm{E}\{{(X_j-{\mathbf{X}}_{-j}^T{\mathbf{\theta }}_j)}^2\}$. When applying the QPC to variable screening, we must estimate two quantities ${\mathbf{\alpha }}_j^0$ and ${\mathbf{\theta }}_j^0$ in advance. However, for ultrahigh dimensional data, the dimensionality of ${\mathbf{X}}_{-j}$ is $p_n-1$, which can still be much bigger than the sample size n. In this situation, it is difficult to obtain the estimators of ${\mathbf{\alpha }}_j^0$ and ${\mathbf{\theta }}_j^0$. On the other hand, it is usually believed that the useful conditional variables are relatively less. Thus, it is reasonable to consider a small subset of $\{k:k\ne j,k=1\dots ,p_n\}$, denoted by ${\mathcal{S}}_j$ in practice. Here, ${\mathcal{S}}_j$ is said to be conditional set with a size smaller than n and it can be specified as the set of previously selected variables and the variables related to the jth predictor, if there is no prior knowledge on it. As a result, Ref. [17] suggested using the following measure to perform variable screening:

$$\begin{array}{c}\hfill {\mathrm{qpcor}}_{\tau }(Y,X_j|{\mathbf{X}}_{{\mathcal{S}}_j})=\frac{\mathrm{E}\{{\psi }_{\tau }(Y-{\mathbf{X}}_{{\mathcal{S}}_j}^T{\mathbf{\alpha }}_j^0)(X_j-{\mathbf{X}}_{{\mathcal{S}}_j}^T{\mathbf{\theta }}_j^0)\}}{\sqrt{\tau (1-\tau ){\sigma }_j^2}},\end{array}$$

(3)

where ${\sigma }_j^2=\mathrm{var}(X_j-{\mathbf{X}}_{{\mathcal{S}}_j}^T{\mathbf{\theta }}_j^0)$, ${\mathbf{\alpha }}_j^0={\mathrm{argmin}}_{{\mathbf{\alpha }}_j}\mathrm{E}\{{\rho }_{\tau }(Y-{\mathbf{X}}_{{\mathcal{S}}_j}^T{\mathbf{\alpha }}_j)\}$ and ${\mathbf{\theta }}_j^0={\mathrm{argmin}}_{{\mathbf{\theta }}_j}\mathrm{E}\{{(X_j-{\mathbf{X}}_{{\mathcal{S}}_j}^T{\mathbf{\theta }}_j)}^2\}$, in which ${\mathbf{X}}_{{\mathcal{S}}_j}={(X_k,k\in {\mathcal{S}}_j)}^T$.

From the definition, one can see that the QPC is just the QC between Y and $X_j$ after removing the confounding effects of conditional variables ${\mathbf{X}}_{{\mathcal{S}}_j}$. Typically, it is through fitting two parametric regression models: one is to fit a linear quantile regression of Y on ${\mathbf{X}}_{{\mathcal{S}}_j}$, and another is on a multivariate linear regression of $X_j$ on ${\mathbf{X}}_{{\mathcal{S}}_j}$. Afterwards, the QPC computes the QC of two residuals that are obtained from these two regression fittings. However, in real applications, the parametric models used to dispel the confounding effects may not be adequate, especially when a nonlinear dependence structure between the response and the predictions is present, which is quite common in high-dimensional data. This motivates us to consider a more flexible and efficient approach to control the influence of the confounding/conditional variables.

2.2. Proposed Method: NQPC-SIS

Without loss of generality, we assume that the predictors $\{X_j,1\le j\le p\}$ are standardized and the response Y satisfies $\tau$-qauntile centered, i.e., $Q_{\tau ,Y}=0$, which is similar to the treatment where the response is centered by mean. Then, we consider the quantile additive model as

$$\begin{array}{c}\hfill Y=m_1(X_1)+m_2(X_2)+\cdots +m_p(X_p)+\epsilon ,\end{array}$$

where the error term satisfies $P(\epsilon <0|\mathbf{X})=\tau$. This means that the conditional $\tau$-quantile of Y given $\mathbf{X}$ is $Q_{\tau ,Y|\mathbf{X}}=m_1(X_1)+m_2(X_2)+\cdots +m_p(X_p)$. We denote by ${\mathcal{M}}_{*}=\{j:m_j(X_j)\ne 0,1\le j\le p\}$ the active set, which indicates the set of indices associated with the nonzero coefficients in the true model and is often assumed to be sparse.

Let $|{\mathcal{S}}_j|$ be the cardinality of a set ${\mathcal{S}}_j$, and $m_{jk}$ and $g_{jk}$, $k\in {\mathcal{S}}_j,$ be ${\ell }_2$-smoothing functions satisfying some conditions. For the identification, we require that $\int m_{jk}(x)\mathrm{d}x=0$ and $\mathrm{E}\{g_{jk}(X_k)\}=0$ for all $j,k$. Set $m_j({\mathbf{X}}_{{\mathcal{S}}_j})={\sum }_{k\in {\mathcal{S}}_j}{m}_{jk}(X_k)$ and $g_j({\mathbf{X}}_{{\mathcal{S}}_j})={\sum }_{k\in {\mathcal{S}}_j}{g}_{jk}(X_k)$. A nonparametric version of QPC (denoted as NQPC) is formulated as

$$\begin{array}{c}\hfill {\varrho }_{\tau }(Y,X_j|{\mathbf{X}}_{{\mathcal{S}}_j})=\frac{\mathrm{E}\{{\psi }_{\tau }(Y-m_j^0({\mathbf{X}}_{{\mathcal{S}}_j}))(X_j-g_j^0({\mathbf{X}}_{{\mathcal{S}}_j}))\}}{\sqrt{\tau (1-\tau ){\sigma }_{j,0}^2}},\end{array}$$

(4)

where ${\sigma }_{j,0}^2=\mathrm{var}(X_j-g_j^0({\mathbf{X}}_{{\mathcal{S}}_j}))$, $m_j^0={\mathrm{argmin}}_{m_j}\mathrm{E}\{{\rho }_{\tau }(Y-m_j({\mathbf{X}}_{{\mathcal{S}}_j}))\}$ and $g_j^0={\mathrm{argmin}}_{g_j}\mathrm{E}\{{(X_j-g_j({\mathbf{X}}_{{\mathcal{S}}_j}))}^2\}$. Suppose we have a dataset: $\{(Y_i,{\mathbf{X}}_i),i=1,\cdots ,n\}$ consisting of n independent copies of $(Y,\mathbf{X})$, where the dimensionality of ${\mathbf{X}}_i$ is $p_n$. Let ${\mathbf{X}}_{i,{\mathcal{S}}_j}$ be the sub-vector of ${\mathbf{X}}_i$ indexed by ${\mathcal{S}}_j$. Then, a sample estimate for NQPC can be given as

$$\begin{array}{c}\hfill {\tilde{\varrho }}_{\tau }(Y,X_j|{\mathbf{X}}_{{\mathcal{S}}_j})=\frac{n^{-1}{\sum }_{i=1}^n{\psi }_{\tau }(Y_i-{\tilde{m}}_j({\mathbf{X}}_{i,{\mathcal{S}}_j}))(X_{ij}-{\tilde{g}}_j({\mathbf{X}}_{i,{\mathcal{S}}_j}))}{\sqrt{\tau (1-\tau ){\tilde{\sigma }}_j^2}},\end{array}$$

(5)

where ${\tilde{\sigma }}_j^2=n^{-1}{\sum }_{i=1}^n{(X_{ij}-{\tilde{g}}_j({\mathbf{X}}_{i,{\mathcal{S}}_j}))}^2$, ${\tilde{m}}_j={\mathrm{argmin}}_{m_j}\frac{1}{n}{\sum }_{i=1}^n{\rho }_{\tau }(Y_i-m_j({\mathbf{X}}_{i,{\mathcal{S}}_j}))$ and ${\tilde{g}}_j={\mathrm{argmin}}_{g_j}\frac{1}{n}{\sum }_{i=1}^n{(X_{ij}-g_j({\mathbf{X}}_{i,{\mathcal{S}}_j}))}^2$. Since $m_{jk}$ and $g_{jk}$ are unknown nonparametric functions, so ${\tilde{m}}_j$ and ${\tilde{g}}_j$ cannot be used, rendering ${\tilde{\varrho }}_{\tau }(Y,X_j|{\mathbf{X}}_{{\mathcal{S}}_j})$ inapplicable. In what follows, we estimate each of $m_{jk}$s and $g_{jk}$s by making use of nonparametric B-spline approximation.

To proceed, we denote $\{B_k(·),k=1,\cdots ,L_n\}$ with $\parallel B_k{\parallel }_{\infty }\le 1$ by a sequence of normalized and centered B-spline basis functions, where $L_n$ is the number of basis functions. Then, according to the theory of B-spline approximation ([29]), for a generic smoothing function m, there exists a vector $\gamma \in {\mathrm{R}}^{L_n}$ such that $m(x)\approx \mathbf{B}{(x)}^T\gamma$, where $\mathbf{B}(·)={(B_1(·),\cdots ,B_{L_n}(·))}^T$. Therefore, there exist vectors ${\mathbf{\alpha }}_{jk}\in {\mathrm{R}}^{L_n}$ and ${\mathbf{\theta }}_{jk}\in {\mathrm{R}}^{L_n}$ such that $m_{jk}(X_k)\approx \mathbf{B}{(X_k)}^T{\mathbf{\alpha }}_{jk}$ and $g_{jk}(X_k)\approx \mathbf{B}{(X_k)}^T{\mathbf{\theta }}_{jk}$. Since $\int m_{jk}(x)\mathrm{d}x=0$ and $\mathrm{E}\{g_{jk}(X_k)\}=0$, it naturally implies that $\mathrm{E}\{\mathbf{B}(X_k)\}=0$ for $k\in {\mathcal{S}}_j$. Write ${\mathbf{\alpha }}_j={(\{{\mathbf{\alpha }}_{jk}^T,k\in {\mathcal{S}}_j\})}^T$, ${\mathbf{\theta }}_j={(\{{\mathbf{\theta }}_{jk}^T,k\in {\mathcal{S}}_j\})}^T$ and ${\mathbf{B}}_j={(\{\mathbf{B}{(X_k)}^T,k\in {\mathcal{S}}_j\})}^T$. Denote by ${\widehat{m}}_j({\mathbf{X}}_{i,{\mathcal{S}}_j})={\mathbf{B}}_{ij}^T{\widehat{\mathbf{\alpha }}}_j$, ${\widehat{g}}_j({\mathbf{X}}_{i,{\mathcal{S}}_j})={\mathbf{B}}_{ij}^T{\widehat{\mathbf{\theta }}}_j$ and ${\widehat{\sigma }}_j^2=n^{-1}{\sum }_{i=1}^n{(X_{ij}-{\widehat{g}}_j({\mathbf{X}}_{i,{\mathcal{S}}_j}))}^2$, where

$$\begin{array}{c}\hfill {\widehat{\mathbf{\alpha }}}_j={\mathrm{argmin}}_{{\mathbf{\alpha }}_j}\frac{1}{n}\sum _{i=1}^n{\rho }_{\tau }(Y_i-{\mathbf{B}}_{ij}^T{\mathbf{\alpha }}_j)\end{array}$$

and

$$\begin{array}{c}\hfill {\widehat{\mathbf{\theta }}}_j={\mathrm{argmin}}_{{\mathbf{\theta }}_j}\frac{1}{n}\sum _{i=1}^n{(X_{ij}-{\mathbf{B}}_{ij}^T{\mathbf{\theta }}_j)}^2,\end{array}$$

where ${\mathbf{B}}_{ij}$ indicates ${\mathbf{B}}_j$ within $\mathbf{B}(X_k)$ being replaced by $\mathbf{B}(X_{ik})$ for $i=1,\cdots ,n$ and $k\in {\mathcal{S}}_j$. Then, it follows that a feasible sample estimate for NQPC is given by

$$\begin{array}{c}\hfill {\widehat{\varrho }}_{\tau }(Y,X_j|{\mathbf{X}}_{{\mathcal{S}}_j})=\frac{n^{-1}{\sum }_{i=1}^n{\psi }_{\tau }(Y_i-{\widehat{m}}_j({\mathbf{X}}_{i,{\mathcal{S}}_j}))(X_{ij}-{\widehat{g}}_j({\mathbf{X}}_{i,{\mathcal{S}}_j}))}{\sqrt{\tau (1-\tau ){\widehat{\sigma }}_j^2}}.\end{array}$$

(6)

Next, we employ the above NQPC estimator as a screening utility in variable screening. To this end, we denote ${\widehat{\mathcal{M}}}_{{\nu }_n}$ to be the selected active set via the screening procedure such that the maximal absolute sample NQPC of the selected variables in ${\widehat{\mathcal{M}}}_{{\nu }_n}$ are greater than a user-specified threshold value ${\nu }_n$. In other words, we can select an active set of variables by

$$\begin{array}{c}\hfill {\widehat{\mathcal{M}}}_{{\nu }_n}=\{j:|\widehat{{\varrho }_{\tau }}(Y,X_j|{\mathbf{X}}_{{\mathcal{S}}_j})|\ge {\nu }_n\mathrm{for}1\le j\le p\}.\end{array}$$

(7)

We name this procedure as the NQPC-based variable screening, abbreviated as NQPC-SIS. In the next section, we will provide some theoretical justification for this approach.

3. Theoretical Properties

To state our theoretical results, we first make some notations. Let $r_n={max}_{1\le j\le p}|{\mathcal{S}}_j|$. Throughout the rest of the paper, for any matrix $\mathbf{A}$, we use $\parallel \mathbf{A}\parallel =\sqrt{{\lambda }_{max}({\mathbf{A}}^T\mathbf{A})}$, ${\parallel \mathbf{A}\parallel }_{\infty }={max}_{i,j}|A_{ij}|$, and ${\lambda }_{min}(\mathbf{A})$ and ${\lambda }_{max}(\mathbf{A})$ to stand for the operator norm, the infinity norm as well as the minimum and maximum eigenvalues for a symmetric matrix $\mathbf{A}$, respectively. In addition, for any vector $\mathbf{a}$, $\parallel \mathbf{a}\parallel =\sqrt{{\sum }_i{a}_i^2}$ means the Euclidean norm.

Denote $u_j=|{\varrho }_{\tau }(Y,X_j|{\mathbf{X}}_{{\mathcal{S}}_j})|$ and ${\widehat{u}}_j=|{\widehat{\varrho }}_{\tau }(Y,X_j|{\mathbf{X}}_{{\mathcal{S}}_j})|$, where ${\varrho }_{\tau }(Y,X_j|{\mathbf{X}}_{{\mathcal{S}}_j})$ is given in Equation (4) and ${\widehat{\varrho }}_{\tau }(Y,X_j|{\mathbf{X}}_{{\mathcal{S}}_j})$ is given in Equation (7). Further, we also denote $u_j^{*}=|{\varrho }_{\tau }^{*}(Y,X_j|{\mathbf{X}}_{{\mathcal{S}}_j})|$, where

$$\begin{array}{c}\hfill {\varrho }_{\tau }^{*}(Y,X_j|{\mathbf{X}}_{{\mathcal{S}}_j})=\frac{\mathrm{E}\{{\psi }_{\tau }(Y-{\mathbf{B}}_j^T{\mathbf{\alpha }}_j^0)(X_j-{\mathbf{B}}_j^T{\mathbf{\theta }}_j^0)\}}{\sqrt{\tau (1-\tau ){\sigma }_j^2}},\end{array}$$

(8)

where ${\sigma }_j^2=\mathrm{var}(X_j-{\mathbf{B}}_j^T{\mathbf{\theta }}_j^0)$, ${\mathbf{\alpha }}_j^0={\mathrm{argmin}}_{{\mathbf{\alpha }}_j}\mathrm{E}\{{\rho }_{\tau }(Y-{\mathbf{B}}_j^T{\mathbf{\alpha }}_j)\}$ and ${\mathbf{\theta }}_j^0={\mathrm{argmin}}_{{\mathbf{\theta }}_j}\mathrm{E}\{{(X_j-{\mathbf{B}}_j^T{\mathbf{\theta }}_j)}^2\}$. Before we establish the uniform convergence of ${\widehat{u}}_j$ to $u_j$, we first investigate the bound of the gap between $u_j$ and $u_j^{*}$, which is helpful to understand the marginal signal level after applying B-spline approximation to the population utility. We need the following conditions:

(B1) 

We assume that $\mathrm{E}\{X_j|{\mathbf{X}}_{{\mathcal{S}}_j}\}=g_j^0({\mathbf{X}}_{{\mathcal{S}}_j})={\sum }_{k\in {\mathcal{S}}_j}{g}_{jk}^0(X_k)$ and ${\mathcal{X}}_k$ denotes the support of covariate $X_k$. There exist some positive constants $C_g$ and $C_m$ such that for any $k\in {\mathcal{S}}_j$,

$$\begin{array}{c}\hfill \underset{1\le j\le p}{max}\underset{x\in {\mathcal{X}}_k}{sup}|g_{jk}^0(x)-{\mathbf{B}}_j{(x)}^T{\mathbf{\theta }}_{jk}^0|\le C_g{L}_n^{-d},\\ \hfill \underset{1\le j\le p}{max}\underset{x\in {\mathcal{X}}_k}{sup}|m_{jk}^0(x)-{\mathbf{B}}_j{(x)}^T{\mathbf{\alpha }}_{jk}^0|\le C_m{L}_n^{-d},\end{array}$$

where d is defined in condition (C1) below.

(B2) 

There exist some positive constants $c_{\sigma ,min},c_{\sigma ,max},{\tilde{c}}_{\sigma ,min},{\tilde{c}}_{\sigma ,max}$ such that

$$\begin{array}{c}\hfill 0<c_{\sigma ,min}\le \underset{1\le j\le p}{max}{\sigma }_j^2\le c_{\sigma ,max}<\infty ,\\ \hfill 0<{\tilde{c}}_{\sigma ,min}\le \underset{1\le j\le p}{max}{\sigma }_{j,0}^2\le {\tilde{c}}_{\sigma ,max}<\infty ,\end{array}$$

where ${\sigma }_j^2$ and ${\sigma }_{j,0}^2$ are given in (4) and (8), respectively.

(B3) 

In a neighborhood of ${\mathbf{B}}_j^T{\mathbf{\alpha }}_j^0$, the conditional density of Y given $(X_j,{\mathbf{X}}_{{\mathcal{S}}_j})$, $f_{Y|(X_j,{\mathbf{X}}_{{\mathcal{S}}_j})}(y)$, is bounded on the support of $(X_j,{\mathbf{X}}_{{\mathcal{S}}_j})$ and uniformly in j.

(B4) 

${min}_{j\in {\mathcal{M}}_{*}}{u}_j\ge C_0{r}_n{n}^{-\kappa }$ for some $C_0>0$ and $0<\kappa <1/2$.

Condition (B1) is imposed on the approximation error condition for nonparametric function in B-spline smoothing literature (e.g., [11,30,31]). Condition (B2) requires variances ${\sigma }_j^2$ and ${\sigma }_{j,0}^2$ to be uniformly bounded. Condition (B3) implies that there exists a finite constant ${\overline{c}}_f>0$ such that for a small $ϵ>0$, ${sup}_{|y-{\mathbf{B}}_j^T{\mathbf{\alpha }}_j^0|<ϵ}{f}_{Y|(X_j,{\mathbf{X}}_{{\mathcal{S}}_j})}(y)\le {\overline{c}}_f$ holds uniformly. Condition (B4) guarantees that the marginal signal of active components in model ${\mathcal{M}}_{*}$ does not vanish. These conditions are similar to those in [17].

Proposition 1.

Under conditions (B1)–(B3), there exists a positive constant $M_{1*}$ such that

$$\begin{array}{c}\hfill u_j-u_j^{*}\le M_{1*}{r}_n{L}_n^{-d},\end{array}$$

In addition, if condition (B4) further holds, then

$$\begin{array}{c}\hfill \underset{j\in {\mathcal{M}}_{*}}{min}{u}_j^{*}\ge C_0\xi r_n{n}^{-\kappa },\end{array}$$

provided that $L_n^{-d}\le C_0(1-\xi )n^{-\kappa }/M_1$ for some $\xi \in (0,1)$.

To establish the sure screening property, we make the following assumptions:

(C1) 

$\{m_{kj}\}$ and $\{g_{kj}\}$ belong to a class of functions $\mathcal{F}$, whose rth derivatives $m_{kj}^{(r)}$ and $g_{kj}^{(r)}$ exist and are Lipschitz of order $\mathbf{\alpha }$,

$$\mathcal{F}=\{b(·):|b^{(r)}(s)-b^{(r)}(t)|\le K|s-t{|}^{\mathbf{\alpha }}\},\mathrm{for}s,t\in [a,b]$$

for some positive constant K, where $[a,b]$ is the support of $X_k$, r is a non-negative integer and $\mathbf{\alpha }\in (0,1]$ such that $d=r+\mathbf{\alpha }>0.5$.

(C2) 

The joint density of $\mathbf{X}$, $f_{\mathbf{X}}$ is bounded by two positive numbers $b_{1f}$ and $b_{2f}$ satisfying $b_{1f}\le f_{\mathbf{X}}\le b_{2f}$. The density of $X_j$, $f_{X_j}$ is bounded away from zero and infinity uniformly in j, that is, there exist two positive constants $c_{1f}$ and $c_{2f}$ such that $c_{1f}\le f_{X_j}(x)\le c_{2f}$.

(C3) 

There exist two positive constants $K_1$ and $K_2$, such that $P(X_j>x|{\mathbf{X}}_{-j})\le K_{1}exp(-K_2^{-1}x)$ for every j.

(C4) 

The conditional density of Y given $\mathbf{X}=\mathbf{x}$, $f_{Y|\mathbf{X}=\mathbf{x}}(y)$, satisfies the Lipschitz condition of first order and $c_{3f}\le f_{Y|\mathbf{X}=\mathbf{x}}(y)\le c_{4f}$ for some positive constants $c_{3f}$ and $c_{4f}$ for any y in a neighborhood of ${\mathbf{B}}_j^T{\mathbf{\alpha }}_j^0$ for $1\le j\le p$.

(C5) 

There exist some positive constants $M_1$ and $M_2$ such that ${sup}_{i,j}|{\mathbf{B}}_{ij}^T{\mathbf{\alpha }}_j^0|\le M_1<\infty$, ${sup}_{i,j}|{\mathbf{B}}_{ij}^T{\mathbf{\theta }}_j^0|\le M_2<\infty$. Furthermore, assume that ${min}_{1\le j\le p}{\sigma }_j^2\ge M_3>0$ for some constant $M_3$.

(C6) 

There exists some constant $\xi \in (0,1)$ such that $L_n^{-d}\le C_0(1-\xi )n^{-\kappa }/M_{1*}$.

Condition (C1) is a smoothness assumption on $\{m_{kj}\}$ and $\{g_{kj}\}$ in nonparametric B-spline-related literature ([7,32]). Condition (C3) is a moment constraint on each of the predictors. Conditions (C2), (C4) and (C5) are similar to those imposed in [17]. Condition (C6) is assumed to ensure the marginal signal level of truly active variables not too weak after B-spline approximation. The above conditions are standard in variable screening literature (e.g., [17,28]).

According to the properties of normalized B-splines and under the conditions (C1) and (C2) (c.f., [33,34]), we can obtain the fact that for each $j=1,\cdots ,p$ and $k=1,\dots ,L_n$, there exist positive constants $C_1,C_2$ and $C_3$ independent of $j,k$ such that

$$\begin{array}{c}\hfill C_1{L}_n^{-1}\le {\lambda }_{min}(\mathrm{E}\{\mathbf{B}(X_j)\mathbf{B}{(X_j)}^T\})\le {\lambda }_{max}(\mathrm{E}\{\mathbf{B}(X_j)\mathbf{B}{(X_j)}^T\})\le C_2{L}_n^{-1},\end{array}$$

(9)

and

$$\begin{array}{c}\hfill \mathrm{E}\{B_k^2(X_j)\}\le C_3{L}_n^{-1}.\end{array}$$

(10)

The following lemma bounds the eigenvalues of the B-spline basis matrix from below and from above. This result extends Lemma 3 of [32] from a fixed dimension to a diverging dimension, which may be crucial to the independent interest of some readers.

Lemma 1.

Suppose that conditions (C1) and (C2) hold, then we have

$$\begin{array}{c}\hfill C_1{\left(\frac{1-{\delta }_0}{2}\right)}^{|{\mathcal{S}}_j|-1}{L}_n^{-1}\le {\lambda }_{min}(\mathrm{E}\{{\mathbf{B}}_j{\mathbf{B}}_j^T\})\le {\lambda }_{max}(\mathrm{E}\{{\mathbf{B}}_j{\mathbf{B}}_j^T\})\le C_2|{\mathcal{S}}_j|L_n^{-1},\end{array}$$

where ${\delta }_0={(1-b_{1f}^2{b}_{2f}^{-2}\zeta )}^{1/2}$ for some constant $0<\zeta <1$.

This result reveals that $r_n$ plays an important role in bounding the eigenvalues of the B-spline basis matrix. When $r_n$ goes to infinity rapidly, the minimum eigenvalue of the basis matrix will degrade to zero very quickly at an exponential rate. However, if the following result holds, then the divergence rate of $r_n$ cannot achieve a polynomial order of n, but can be of an order of $logn$.

Theorem 1.

Suppose that conditions (B1)–(B5) and (C1)–(C5) hold and assume that $a_0^{-2r_n}{L}_n/n^{1-2\kappa }=o(1)$ and $a_0^{-2r_n}{r}_n^3{L}_n{n}^{-\kappa }=o(1)$ are satisfied.

(i) 

For any $C>0$, then there exist some positive constants $c_6^{*},c_{14}^{*}$ such that, for $0<\kappa <1/2$ and sufficiently large n,

$$\begin{array}{ccc}& & P\left(\underset{1\le j\le p_n}{max}|{\widehat{u}}_j-u_j^{*}|\ge C{r}_n{n}^{-\kappa }\right)\hfill \\ & & \le p_n\{7exp(-c_6^{*}{a}_0^{2r_n}{r}_n^2{n}^{1-4\kappa })+[116{(r_n{L}_n)}^2+60r_n{L}_n+10]exp(-c_{14}^{*}{a}_0^{2r_n}{L}_n^{-3}{n}^{1-2\kappa })\},\hfill \end{array}$$

where $a_0=(1-{\delta }_0)/2$ and ${\delta }_0$ is given in Lemma 1.

(ii) 

In addition, if condition (C6) is further satisfied, by choosing ${\nu }_n={\tilde{C}}_0{r}_n{n}^{-\kappa }$ with ${\tilde{C}}_0\le C_0\xi /2$, we have

$$\begin{array}{ccc}\hfill P\left({\mathcal{M}}_{*}\subset {\widehat{\mathcal{M}}}_{{\nu }_n}\right)& & \ge 1-s_n\{7exp(-c_6^{*}{a}_0^{2r_n}{r}_n^2{n}^{1-4\kappa })\hfill \\ & & +[116{(r_n{L}_n)}^2+60r_n{L}_n+10]exp(-c_{14}^{*}{a}_0^{2r_n}{L}_n^{-3}{n}^{1-2\kappa })\}\hfill \end{array}$$

for sufficiently large n, where $s_n=|{\mathcal{M}}_{*}|$.

The above establishes the sure screening property that all the relevant variables can be recruited with probability going to one in the final model. The probability bound in the property is free of $p_n$, but depends on $r_n$ and the number of basis functions $L_n$. Though this ensures that NQPC-SIS retains all important predictors with high probability, the noisy variables can be included by NQPC-SIS. Ideally, this can be realized by the choice of ${\nu }_n$, according to Theorem 1 and by setting ${max}_{j\notin {\mathcal{M}}_{*}}|{\varrho }_{\tau }^{*}(Y,X_j|{\mathbf{X}}_{{\mathcal{S}}_j})|=o(r_n{n}^{-\kappa })$, to achieve the selection consistency, i.e.,

$$\begin{array}{c}\hfill P({\mathcal{M}}_{*}={\widehat{\mathcal{M}}}_{{\nu }_n})\to 1\end{array}$$

when n is sufficiently large. This property can also be achieved by Theorem 1 and by assuming that ${\varrho }_{\tau }^{*}(Y,X_j|{\mathbf{X}}_{{\mathcal{S}}_j})=0$ for $j\notin {\mathcal{M}}_{*}$. However, this would be too restrictive to check in practice. Similar to [17], we may assume that ${\sum }_{j=1}^p{u}_j^{*}=O(n^{\varsigma })$ for some $\varsigma >0$ to control the false selection rate. With this condition, we can obtain the following property to control the size of the selected model.

Theorem 2.

Under the conditions of Theorem 1 and by choosing ${\nu }_n={\tilde{C}}_0{r}_n{n}^{-\kappa }$ with ${\tilde{C}}_0\le C_0\xi /2$ and if ${\sum }_{j=1}^p{u}_j^{*}=O(n^{\varsigma })$ for some $\varsigma >0$, then for some positive constant $C_{*}$, there exist some constants ${\tilde{c}}_6^{*},{\tilde{c}}_{14}^{*}$ such that

$$\begin{array}{ccc}\hfill P(|{\widehat{\mathcal{M}}}_{{\nu }_n}|\le C_{*}{r}_n^{-1}{n}^{\kappa +\varsigma })& & \ge 1-p_n\{7exp(-{\tilde{c}}_6^{*}{a}_0^{2r_n}{r}_n^2{n}^{1-4\kappa })\hfill \\ & & +[116{(r_n{L}_n)}^2+60r_n{L}_n+10]exp(-{\tilde{c}}_{14}^{*}{a}_0^{2r_n}{L}_n^{-3}{n}^{1-2\kappa })\}\hfill \end{array}$$

for sufficiently large n.

This theorem reveals that after an application of the NQPC-SIS, the dimensionality can be reduced from an exponential order to a polynomial size of n at the same time retaining all the important predictors with probability approaching one.

4. Algorithm for NQPC-SIS

To make the NQPS-SIS practically applicable, for each $X_J$, we need to specify the conditional set ${\mathcal{S}}_j$. We note that a sequential test was developed to identify ${\mathcal{S}}_j$ in [17] via an application of the Fisher’s Z-transformation [35] and partial correlation. In this section, we provide a two-stage procedure based on nonparametric additive quantile regression model, which can be viewed as a complementary to [17].

To reduce the computational burden, we first apply the quantile-adaptive model-free feature screening (Qa-SIS) proposed by [13] to select a subset from $\{X_j,1\le j\le p_n\}$, denoted by ${\widehat{\mathcal{M}}}_{\text{Qa-SIS}}$ with $|{\widehat{\mathcal{M}}}_{\text{Qa-SIS}}|=\lfloor 0.5n{L}_n^{-1}/log(n{L}_n^{-1})\rfloor +1$, where $L_n$ is the number of basis functions used in Qa-SIS and $\lfloor a\rfloor$ denotes the largest integer not exceeding a. Second, for each $X_j$, if $X_j\in {\widehat{\mathcal{M}}}_{\text{Qa-SIS}}$, we set ${\mathcal{C}}_j=\{X_k|X_k\in {\widehat{\mathcal{M}}}_{\text{Qa-SIS}},k\ne j\}$, otherwise ${\mathcal{C}}_j=\{X_k|X_k\in {\widehat{\mathcal{M}}}_{\text{Qa-SIS}},k\ne |{\widehat{\mathcal{M}}}_{\text{Qa-SIS}}|\}$. Thus, $|{\mathcal{C}}_j|=\lfloor 0.5n{L}_n^{-1}/log(n{L}_n^{-1})\rfloor$. Third, we carry out a variable selection with SCAD penalty [2] based on additive quantile regression model for data set $\{(X_{ij},X_{i{\mathcal{C}}_j}),i=1,\cdots ,n\}$ and then a small reduced subset is obtained, denoted by ${\mathcal{C}}_j^v$. Such a two-stage procedure can help to find the conditional subset for the jth variable and will be incorporated in the following algorithm. With a slight abuse of notation, we use $d_n$ to denote the screening threshold parameter of the NQPC-SIS, in other words, for the NQPC-SIS, we select $d_n$ covariates that correspond to the first $d_n$ largest NQPCs.

Algorithm 1 has the same spirit as the QPCS algorithm of [17], who demonstrated empirically that the QPCS algorithm outperforms their QTCS and QFR algorithms. In the implementation, we choose $d_n^{*}=\lfloor 0.5n{L}_n^{-1}/log(n{L}_n^{-1})\rfloor$ and $d_n=\lfloor n/logn\rfloor$, which does not exclude other choice. According to our limited simulation experience, this choice works satisfactorily. The values of $d_n^{*}$ and $r_n$ we take on cannot be too large, due to the use of B-spline basis approximations. Theoretically, we need to specify $d_n^{*}$ such that $d_n^{*}\le r_n$, while it is sufficient to require $L_n{d}_n^{*}<n$ practically.

Algorithm 1 The implementation of NQPC-SIS.

1: Given $d_n$, we set a pre-specified number $d_n^{*}\le d_n$ and an initial set ${\mathcal{A}}^{(0)}=\varnothing$. 2: For $k=1,\dots ,d_n^{*}$, (2a) update ${\mathcal{S}}_j={\mathcal{A}}^{(k-1)}\cup {\mathcal{C}}_j^v$; (2b) update ${\mathcal{A}}^{(k)}={\mathcal{A}}^{(k-1)}\cup \{j^{*}\}$, where the variable index $j^{*}$ is defined by $$\begin{array}{c}\hfill j^{*}={\mathrm{argmax}}_{j\notin {\mathcal{A}}^{(k-1)}}|{\widehat{\rho }}_{\tau }(Y,X_j|X_{{\mathcal{S}}_j})|.\end{array}$$ 3: For $k=d_n^{*}+1,\dots ,d_n$, (3a) update ${\mathcal{S}}_j={\mathcal{A}}^{(d_n^{*})}\cup {\mathcal{C}}_j^v$; (3b) update ${\mathcal{A}}^{(k)}={\mathcal{A}}^{(k-1)}\cup \{j^{*}\}$, where the variable index $j^{*}$ is such that $$\begin{array}{c}\hfill j^{*}={\mathrm{argmax}}_{j\notin {\mathcal{A}}^{(k-1)}}|{\widehat{\rho }}_{\tau }(Y,X_j|X_{{\mathcal{S}}_j})|.\end{array}$$ 4: Repeat Step 3 until $k\ge d_n$. The final selected set is denoted as $\widehat{\mathcal{M}}$.

5. Numerical Studies

5.1. Simulations

In this subsection, we conduct some simulation studies to examine the finite sample performance of the proposed NQPC-SIS. In order to evaluate the performance, we employ three criteria: the minimum model size (MMS), i.e., the smallest number of covariates that contain all the active variables, its robust standard deviation (RSD), and the proportion of all the active variables selected ($\mathcal{P}$) with the screening threshold parameter being specified as $d_n=\lfloor n/logn\rfloor$. Throughout this subsection, we adopt the following simulation settings: the sample size $n=200$, the number of basis $L_n=\lfloor n^{1/5}\rfloor +1$, and the dimensionality $p_n=1000$. We simulate the random error $\epsilon$ from two distributions: $N(0,1)$ and $t(3)$, respectively. Three quantile levels $\tau =0.2,0.5,0.8$ are considered in all situations. For each simulation scenario, all the results are obtained over $N=200$ replications.

Example 1.

Let $\mathbf{X}={(X_1,\cdots ,X_{p_n})}^T$ be a $p_n$-dimensional random vector having a multivariate normal distribution with mean zero and covariance matrix $\Sigma ={({\sigma }_{jk})}_{1\le j,k\le p_n}$, where ${\sigma }_{jj}=1$ and ${\sigma }_{j,k}=\rho ,j\ne k$ except that ${\sigma }_{i4}={\sigma }_{4j}=\sqrt{\rho }$. Generate the response as:

$$\begin{array}{c}\hfill Y=\beta X_1+\beta X_2+\beta X_3-3\beta \sqrt{\rho }{X}_4+\epsilon .\end{array}$$

It is easily observed that the marginal Pearson’s correlation between $X_4$ and Y is zero. We take $\rho =0.5,0.8$ and set $\beta =2.5(1+|\tau -0.5|)$ to incorporate the quantile information.

Example 2.

We follow the simulation model of [17] and generate the response as

$$\begin{array}{c}\hfill Y=\beta X_1+\beta X_2+\beta X_3-3\beta \sqrt{\rho }{X}_4-0.25\beta X_5+\epsilon ,\end{array}$$

where $\beta ,\rho ,$ and $\mathbf{X}$ are defined as in Example 1 except that ${\sigma }_{i5}={\sigma }_{5j}=0$ such that $X_5$ is uncorrelated with $X_j,j\ne 5$.

Example 3.

We simulate the response from the following nonlinear model:

$$\begin{array}{c}\hfill Y=3g_1(X_1)+3g_2(X_2)+3g_3(X_3)+3g_4(X_4)+3g_5(X_5)+\epsilon ,\end{array}$$

where $g_1(x)=1.5x,g_2(x)=2x(2x-1),g_3(x)=sin(2\pi x)/(2sin(2\pi x)),g_4(x)=sin(2\pi x),g_5(x)=e^{x-0.5}$. The covariates $\mathbf{X}=(X_1,\cdots ,X_{p_n})$ are simulated from a random-effects model $X_j=\frac{W_j+tU}{1+t},j=1,\cdots ,p_n$, where $W_j$s and U are iid $Unif(0,1)$. We consider two cases of $t=1$ and $t=2$, corresponding to $\mathrm{corr}(X_j,X_k)=0.5$ and $0.8$ for $j\ne k$, respectively.

Example 4.

We consider the same model as that in Example 3, with exception that $X_2$ and $X_5$ are replaced by $X_2=cos(2\pi X_6)+ϵ$ and $X_5={(X_1-0.5)}^2+ϵ$, where $ϵ\sim N(0,1)$ and is independent of ε, the error in the model in Example 3.

The simulation results of Examples 1–4 are shown in

Table 1. Simulation results for Example 1 when n = 200.

				$\mathit{\tau }=0.2$		$\mathit{\tau }=0.5$		$\mathit{\tau }=0.8$
$\mathit{\epsilon }$	$\mathit{\rho }$	Method	${\mathit{s}}_{\mathit{n}}$	MMS(RSD)	$\mathcal{P}$	MMS(RSD)	$\mathcal{P}$	MMS(RSD)	$\mathcal{P}$
$N(0,1)$	0.5	SIS	4	455.5(319.3)	0	437(330.3)	0	434.5(372)	0
		NIS	4	451(456.5)	0.025	506(421.3)	0	486.5(390.5)	0
		Qa-SIS	4	466(392.5)	0.02	466.5(375.5)	0.01	490.5(382.3)	0.01
		QPC-SIS	4	4(0)	1	4(0)	1	4(0)	1
		NQPC-SIS	4	4(0)	0.995	4(0)	1	4(0)	1
	0.8	SIS	4	444.5(141.3)	0	458(161.8)	0	452.5(188)	0
		NIS	4	489.5(274.5)	0	518.5(274)	0	511(285.8)	0
		Qa-SIS	4	522(372.3)	0.01	510.5(358)	0	560.5(292.8)	0
		QPC-SIS	4	5(2)	0.99	4(1)	1	5(2)	0.98
		NQPC-SIS	4	6(3)	0.96	4(2)	0.99	6(3)	0.96
$t(3)$	0.5	SIS	4	434.5(352.8)	0.005	475(343.3)	0	472.5(366)	0
		NIS	4	492.5(347.5)	0.01	501.5(415)	0	555.5(352.3)	0
		Qa-SIS	4	510.5(390.3)	0.005	481(463.3)	0.015	541.5(460.3)	0.01
		QPC-SIS	4	4(0)	1	4(0)	1	4(0)	1
		NQPC-SIS	4	4(0)	0.995	4(0)	1	4(0)	1
	0.8	SIS	4	453(135.8)	0	468(200.5)	0	473(283.5)	0
		NIS	4	535.5(288.3)	0	507(253.3)	0	507.5(368.3)	0
		Qa-SIS	4	597.5(329.3)	0	578.5(374)	0	591.5(366.8)	0.005
		QPC-SIS	4	6(3)	0.915	5(2)	0.975	6(2)	0.945
		NQPC-SIS	4	6.5(6)	0.84	5(3)	0.955	6(5.3)	0.855

,

Table 2. Simulation results for Example 2 when $n=200$.

				$\mathit{\tau }=0.2$		$\mathit{\tau }=0.5$		$\mathit{\tau }=0.8$
$\mathit{\epsilon }$	$\mathit{\rho }$	Method	${\mathit{s}}_{\mathit{n}}$	MMS(RSD)	$\mathcal{P}$	MMS(RSD)	$\mathcal{P}$	MMS(RSD)	$\mathcal{P}$
$N(0,1)$	0.5	SIS	5	439.5(359.5)	0	477(319)	0	427(324.8)	0
		NIS	5	522(362)	0.005	566(429.8)	0	507.5(392)	0
		Qa-SIS	5	542.5(400.3)	0	565.5(351.5)	0	554(340.3)	0
		QPC-SIS	5	5(0)	1	5(0)	1	5(0)	1
		NQPC-SIS	5	5(0)	1	5(0)	1	5(0)	1
	0.8	SIS	5	436(111.3)	0	479.5(232.8)	0	466.5(219.3)	0
		NIS	5	523.5(246.8)	0	556.5(265.8)	0	527(286.3)	0
		Qa-SIS	5	557.5(376.5)	0	604(358.8)	0	542(363.8)	0
		QPC-SIS	5	6(2)	0.97	6(2)	0.97	7(3)	0.93
		NQPC-SIS	5	7(2)	0.9	6(2)	0.945	7(3)	0.9
$t(3)$	0.5	SIS	5	478.5(347)	0	451.5(308.8)	0	483.5(322.8)	0
		NIS	5	535.5(384.3)	0	545.5(317.8)	0.005	508(353)	0
		Qa-SIS	5	597.5(389.5)	0.005	568.5(341.5)	0	593.5(435.8)	0.005
		QPC-SIS	5	5(0)	0.99	5(0)	1	5(0)	0.995
		NQPC-SIS	5	5(1)	0.985	5(0)	0.995	5(1)	0.975
	0.8	SIS	5	468(286.3)	0	477(238.5)	0	466(229)	0
		NIS	5	530.5(324.8)	0	532.5(321.3)	0	525.5(245.5)	0
		Qa-SIS	5	655(391.8)	0	590.5(361.5)	0	591.5(374.8)	0
		QPC-SIS	5	7(21.3)	0.765	7(3)	0.88	8(44.8)	0.74
		NQPC-SIS	5	8(81.8)	0.63	7(7.3)	0.82	11(136.5)	0.57

,

Table 3. Simulation results for Example 3 when $n=200$.

				$\mathit{t}=1$		$\mathit{t}=2$
$\mathit{\epsilon }$	$\mathit{\tau }$	Method	${\mathit{s}}_{\mathit{n}}$	MMS(RSD)	$\mathcal{P}$	MMS(RSD)	$\mathcal{P}$
$N(0,1)$	$\tau =0.2$	Qa-SIS	5	695(461.5)	0.005	799.5(326.5)	0
		QPC-SIS	5	46(165.3)	0.485	298.5(477.5)	0.075
		NQPC-SIS	5	7(13)	0.82	160(345.8)	0.25
	$\tau =0.5$	SIS	5	492(805)	0.155	1000(1)	0
		NIS	5	533(608)	0.035	798(391)	0.005
		Qa-SIS	5	761.5(460)	0.005	762.5(353)	0
		QPC-SIS	5	8.5(33)	0.75	211(372.8)	0.145
		NQPC-SIS	5	5(0)	0.96	29(103.8)	0.56
	$\tau =0.8$	Qa-SIS	5	599(508)	0.01	708(348.3)	0
		QPC-SIS	5	47(213)	0.46	393(400)	0.015
		NQPC-SIS	5	6(7)	0.86	156.5(388)	0.23
$t(3)$	$\tau =0.2$	Qa-SIS	5	689.5(434)	0	794(345)	0
		QPC-SIS	5	85.5(239.5)	0.28	548.5(518.5)	0.025
		NQPC-SIS	5	46.5(181)	0.46	487(551)	0.055
	$\tau =0.5$	SIS	5	626.5(767.8)	0.1	999(6)	0
		NIS	5	560.5(574.5)	0.025	742(430.3)	0
		Qa-SIS	5	673.5(449.3)	0.005	751.5(358.5)	0
		QPC-SIS	5	21.5(89.3)	0.56	331.5(467.8)	0.06
		NQPC-SIS	5	6(5)	0.855	136.5(388)	0.21
	$\tau =0.8$	Qa-SIS	5	583(458.5)	0.015	711.5(382.3)	0
		QPC-SIS	5	108.5(303)	0.3	623(448)	0.005
		NQPC-SIS	5	28.5(136.3)	0.52	413(489.5)	0.03

and

Table 4. Simulation results for Example 4 when $n=200$.

				$\mathit{t}=1$		$\mathit{t}=2$
$\mathit{\epsilon }$	$\mathit{\tau }$	Method	${\mathit{s}}_{\mathit{n}}$	MMS(RSD)	$\mathcal{P}$	MMS(RSD)	$\mathcal{P}$
$N(0,1)$	$\tau =0.2$	Qa-SIS	5	735(448.5)	0.005	681.5(388.5)	0
		QPC-SIS	5	647.5(425)	0.005	746.5(329.3)	0
		NQPC-SIS	5	5(1)	0.945	86(334)	0.385
	$\tau =0.5$	SIS	5	793(268.3)	0	846.5(252.8)	0
		NIS	5	765.5(298.8)	0	896.5(225.8)	0
		Qa-SIS	5	749.5(274.8)	0	818(301.5)	0
		QPC-SIS	5	717.5(326.8)	0	805.5(254.3)	0
		NQPC-SIS	5	5(0)	1	8(60.3)	0.705
	$\tau =0.8$	Qa-SIS	5	836(274)	0	867.5(243.8)	0
		QPC-SIS	5	798.5(248.3)	0	811.5(249.3)	0
		NQPC-SIS	5	5(1)	0.985	61(355.8)	0.44
$t(3)$	$\tau =0.2$	Qa-SIS	5	716.5(374.3)	0	703(375.3)	0
		QPC-SIS	5	603(422)	0.01	743.5(295.3)	0
		NQPC-SIS	5	7(22.5)	0.78	317(592.3)	0.15
	$\tau =0.5$	SIS	5	786.5(261.8)	0	869.5(301.3)	0
		NIS	5	779(285.5)	0	833.5(259.8)	0
		Qa-SIS	5	754.5(324.5)	0	800(255.8)	0
		QPC-SIS	5	755.5(379.5)	0	825(296.3)	0
		NQPC-SIS	5	5(0)	0.99	61.5(302.3)	0.435
	$\tau =0.8$	Qa-SIS	5	819.5(255.8)	0	869(241.3)	0
		QPC-SIS	5	795(249.3)	0	847(260.3)	0
		NQPC-SIS	5	6(9)	0.835	375(576)	0.14

, respectively. The results in Table 1 show that when the true relation between the response and covariates in the model is linear, the SIS, NIS and Qa-SIS methods fail to work. However, when comparing to those methods, we can see that both QPC-SIS and NQPC-SIS with $\tau =0.5$ work reasonably well, although the QPC-SIS slightly outperform our NQPC-SIS when $\rho =0.8$. This is expected because the QPC works for the model with linear relationship between the covariates. A similar observation can be drawn in Table 2 for Example 2, which is also a linear model, albeit the difference that $X_5$ and $X_j,j\ne 5$ are independent in Example 2. The results in Table 3 indicate that when the relationship between Y and $\mathbf{X}$ is nonlinear and the relationship between covariates is linear, our proposed NQPC-SIS performs best and then followed by QPC-SIS. From Table 4, we can see that when the relationship between Y and $\mathbf{X}$ is nonlinear and there also exists a nonlinear relationship among $\mathbf{X}$, NQPC-SIS works most satisfactorily and is much better than Qa-SIS and QPC-SIS in terms of both MMS and selection rate $\mathcal{P}$.

In addition, the simulation results of QPC-SIS and NQPC-SIS for Examples 1–4 with $\rho =0.9$ and $\tau =0.5$ are reported in

Table 5. Simulation results for Examples 1 to 4 when $\rho =0.9$ and $\tau =0.5$.

		$\mathit{n}=200$				$\mathit{n}=400$
		$\mathbf{\epsilon }\sim \mathbf{N}(\mathbf{0},\mathbf{1})$		$\mathbf{\epsilon }\sim \mathbf{t}(\mathbf{3})$		$\mathbf{\epsilon }\sim \mathbf{N}(\mathbf{0},\mathbf{1})$		$\mathbf{\epsilon }\sim \mathbf{t}(\mathbf{3})$
	Method	MMS(RSD)	$\mathcal{P}$	MMS(RSD)	$\mathcal{P}$	MMS(RSD)	$\mathcal{P}$	MMS(RSD)	$\mathcal{P}$
Example 1	QPC-SIS	6(2)	0.94	9(64)	0.665	5(2)	1	5(2)	0.985
	NQPC-SIS	6(3)	0.88	30.5(174.3)	0.52	5(2)	0.99	6(2)	0.93
Example 2	QPC-SIS	7(2)	0.9	25.5(204.5)	0.525	6(2)	1	7(2)	0.985
	NQPC-SIS	8(9.25)	0.825	44.5(178.8)	0.465	6(2)	0.995	7(2)	0.95
Example 3	QPC-SIS	608.5(471.5)	0	687.5(400)	0.005	253(408.3)	0.175	399.5(521)	0.105
	NQPC-SIS	341(505.8)	0.09	527.5(496.5)	0.025	19(75)	0.705	69.5(247.8)	0.475
Example 4	QPC-SIS	770.5(311.5)	0	802.5(232.8)	0	752.5(312)	0	784(258.3)	0
	NQPC-SIS	301.5(484.3)	0.11	433(554.8)	0.075	11(44)	0.76	38.5(142.8)	0.58

. It can be observed from Table 5 that when the sample size increases from 200 to 400, the performance of QPC-SIS and NQPC-SIS are improved by much, although QPC-SIS and NQPC-SIS perform very competitively in Examples 1 and 2, while NQPC-SIS performs significantly better than QPC-SIS in Examples 3 and 4. These evidences indicate the effectiveness and usefulness of our NQPC-SIS.

As suggested by one anonymous reviewer, we add one more simulation to compare our NQPC-SIS with the following two approaches: (a) QC-SIS, which is the screening method based on quantile correlation, but simply ignores the effect of conditional variables on the response, and (b) RFQPC-SIS, which is a procedure very similarly to our NQPC-SIS, yet removes the effect of conditional variables through fitting Random Forest models. We examine the performance of these three approaches under $\tau =0.5$ and $n=200$ for Examples 1 to 4, where RFQPC-SIS is a variant of the NQPC method and implemented with randomForest in R package “randomForest”. Note that RFQPC-SIS requires $2(p-|{\mathcal{A}}^{(k)}|)$ random forest regressions in the k-th iteration, which is highly computationally intensive. Here, we evaluate NQPC-SIS, QC-SIS and RFQPC-SIS using effective model size (EMS) and $\mathcal{P}$, where EMS indicates the average of true variables contained in the first $d_n=\lfloor n/log(n)\rfloor$ variables selected from 200 replicate experiments. The results are reported in

Table 6. Simulation results for Examples 1 to 4 when $n=200$ and $\tau =0.5$.

		$\mathit{\rho }=0.5$				$\mathit{\rho }=0.8$
		$\mathbf{\epsilon }\sim \mathbf{N}(\mathbf{0},\mathbf{1})$		$\mathbf{\epsilon }\sim \mathbf{t}(\mathbf{3})$		$\mathbf{\epsilon }\sim \mathbf{N}(\mathbf{0},\mathbf{1})$		$\mathbf{\epsilon }\sim \mathbf{t}(\mathbf{3})$
	Method	EMS(SD)	$\mathcal{P}$	EMS(RSD)	$\mathcal{P}$	EMS(SD)	$\mathcal{P}$	EMS(SD)	$\mathcal{P}$
Example 1	QC-SIS	2.985(0.122)	0	2.925(0.346)	0	2.570(0.969)	0	2.340(1.162)	0
	RFQPC-SIS	3.995(0.071)	0.995	3.975(0.157)	0.975	3.675(0.470)	0.675	3.460(0.500)	0.46
	NQPC-SIS	4(0)	1	4(0)	1	3.985(0.122)	0.985	3.930(0.256)	0.93
Example 2	QC-SIS	3.565(0.646)	0	3.490(0.657)	0	3.355(1.147)	0	3.255(1.280)	0
	RFQPC-SIS	4.745(0.437)	0.745	4.630(0.504)	0.64	4.475(0.609)	0.535	4.215(0.649)	0.335
	NQPC-SIS	5(0)	1	5(0)	1	4.960(0.221)	0.965	4.785(0.447)	0.8
Example 3	QC-SIS	2.470(0.862)	0.03	2.575(0.894)	0	1.505(0.567)	0	1.460(0.592)	0
	RFQPC-SIS	4.945(0.229)	0.945	4.730(0.788)	0.75	4.275(0.808)	0.465	3.535(0.175)	0.175
	NQPC-SIS	4.955(0.208)	0.955	4.785(0.424)	0.79	4.320(0.825)	0.49	3.715(1.109)	0.265
Example 4	QC-SIS	1.775(0.613)	0	1.790(0.720)	0	1.685(0.639)	0	1.690(0.629)	0
	RFQPC-SIS	3.045(0.739)	0.02	3.040(0.788)	0.045	2.455(0.616)	0	2.350(0.528)	0
	NQPC-SIS	5(0)	1	4.980(0.140)	0.98	4.590(0.731)	0.715	4.205(0.864)	0.465

, showing that our NQPC-SIS still performs the best and is followed by RFQPC-SIS. Moreover, the computational cost of NQPC-SIS is much less than that of RFQPC-SIS.

5.2. An Application to Breast Cancer Data

In this subsection, we apply the proposed NQPC-SIS to breast cancer data with a high lethality rate, which is reported by [36]. The data consists of 19,672 gene expression and 2,149 CGH measurements from 89 cancer patient samples, which is available at https://github.com/bnaras/PMA/blob/master/data/breastdata.rda (accessed on 18 June 2021). Our interest here is to detect the genes that have the most impact on comparative genomic hybridization (CGH) measurements. A similar purpose was achieved in [25,37]. Following [37], we consider the first principal component of 136 CGH measurements as the response Y and the remaining 18,672 gene probes as the explanatory variables $\mathbf{X}$. We implement the two stage procedure for the sake of comparison, where a variable screening method is implemented in the first stage and a predictive regression model is conducted in the second stage. To this end, we select $d_n=\lfloor n/log(n)\rfloor$ variables in the first stage using one of the screening methods: SIS, NIS, Qa-SIS, QPC-SIS and NQPC-SIS, as mentioned in the simulation study. In the second stage, we randomly select 80% sample data as the training set, and the remaining $20\%$ sample as the test set. Then, we apply one machine learning method, regression tree, to the dimension-reduced data to examine the finite sample performance on the test set. We use the command M5P in R package “RWeka” for implementing the regression tree method. We use the mean of absolute prediction error (MAPE), defined as

$$\begin{array}{c}\hfill \mathrm{MAPE}=\frac{1}{n^{(test)}}\sum _{i=1}^{n^{(test)}}\left|Y_i^{(test)}-{\widehat{Y}}_i^{(test)}\right|,\end{array}$$

as our evaluation index, where $n^{(test)}$ is the number of observations in the training set and ${\widehat{Y}}_i^{(test)}$ is the predicted value of Y at the observation $x_i$ in the test set. We repeat the above procedure 500 times and report the mean and standard deviation of 500 MAPEs in

Table 7. Prediction results for the real data on the test set, where the standard deviation is given in the parenthesis.

Method	$\mathit{\tau }=0.3$	$\mathit{\tau }=0.4$	$\mathit{\tau }=0.5$	$\mathit{\tau }=0.6$	$\mathit{\tau }=0.7$
SIS	-	-	0.8202(0.1362)	-	-
NIS	-	-	0.8254(0.1348)	-	-
Qa-SIS	0.8318(0.1472)	0.8261(0.1446)	0.8375(0.1448)	0.8612(0.1541)	0.8431(0.1512)
QPC-SIS	0.7269(0.1267)	0.7989(0.1458)	0.8347(0.1471)	0.6495(0.1155)	1.0240(0.1732)
NQPC-SIS	0.7629(0.1356)	0.6488(0.1232)	0.6742(0.1285)	0.8156(0.1643)	0.7802(0.1315)

. According to the results in Table 7, we can observe that the NQPC-SIS outperforms both the SIS, NIS and Qa-SIS. Typically, our NQPC-SIS produces the lowest prediction error (MAPE) among these methods when $\tau =0.4$, $\tau =0.5$ and $\tau =0.7$. Moreover, we also note that the QPC-SIS performs better than our NQPC-SIS at $\tau =0.3$ and $\tau =0.6$, but worse than our method at other three quantile levels. Qa-SIS performs worst among these methods. This evidence supports that the proposed NPQC-SIS in this paper works well for this real data.

6. Concluding Remarks

In this paper, we proposed a nonparametric quantile partial correlation-based variable screening approach (NQPC-SIS), which can be viewed as an extension of the QPC-SIS proposed in [17] from a parametric framework to the nonparametric situation. Our proposed NQPC-SIS enjoys the sure independence screening property under some mild technical conditions. Furthermore, an algorithm of NQPC-SIS for implementation is provided for users. Extensive numerical experiments including simulations and real-world data analysis are carried out for illustration. The numerical results showed that our NQPC-SIS works fairly well especially when the relationship between variables is highly nonlinear.