Model Selection for High Dimensional Nonparametric Additive Models via Ridge Estimation

Abstract

In ultrahigh dimensional data analysis, to keep computational performance well and good statistical properties still working, nonparametric additive models face increasing challenges. To overcome them, we introduce a methodology of model selection for high dimensional nonparametric additive models. Our approach is to propose a novel group screening procedure via nonparametric smoothing ridge estimation (GRIE) to find the importance of each covariate. It is then combined with the sure screening property of GRIE and the model selection property of extended Bayesian information criteria (EBIC) to select the suitable sub-models in nonparametric additive models. Theoretically, we establish the strong consistency of model selection for the proposed method. Extensive simulations and two real datasets illustrate the outstanding performance of the GRIE-EBIC method.

1. Introduction

With the advances in information technology, high-dimensional data exists in various fields such as biology, chemistry, economics, finance, genetics, neuroscience, etc. A common assumption of sparse is that only a few features are truly related to the response. Following that, plenty of variable selection approaches based on regularized M-estimation have been developed, which include but are not limited to Lasso by [1],  SCAD by [2], Dantzig selector by [3], and  MCP by [4]. However, there always exist two limitations in the above-penalized methods. One is the big burden for computation, and the other is the unstable performance for variable selection in high-dimensional situations [5].

To avoid the mentioned limitations, correlation ranking becomes one of the most popular ways to rapidly reduce the dimensionality of feature space. Fan and Lv [6] proposed the sure independence screening (SIS) by utilizing the marginal Pearson correlation between predictor and response for gaussian linear regression. Fan et al. [7] extended the idea of Pearson correlation ranking to marginal smooth estimation strength ranking and proposed the nonparametric independence screening (NIS) method. Meanwhile, Zhu et al. [8] considered the marginal correlation between the predictor and the conditional cumulative density function of response and developed the model-free screening method. However, in  practice, there exist strong correlations between the predictors, which may lead to important predictors being jointly correlated to the response. Hence, the marginal correlation ranking process may miss some important variables. To decrease the effect of correlation between the predictors, some forward variable screening methods based on the prediction rankings were introduced. Wang [9] ranked the residuals of the predictor and proposed the forward regression (FR) algorithm. Cheng et al. [10] applied the forward regression to high dimensional varying coefficient models and proposed the forward-BIC screening method. Zhong et al. [11] further extended the forward regression to ultrahigh-dimensional nonparametric additive models. Based on the cumulative divergence (CD), Zhou et al. [12] proposed a forward screening procedure that considered the joint effects among covariates in the feature screening process.

Next, let us turn to the specific model. In this paper, we are interested in the model of nonlinear regression. It is well known that if there exists extensive nonlinear independence between response and predictors, traditional (partial) linear models can not detect nonlinear independence. Although the nonlinear regression could capture the nonlinear independence accurately, the nonlinear regression suffers from the curse of dimensionality and heavy computational burden in high dimensions. To model them simplify, here, we consider the nonparametric additive models. The nonparametric additive models were introduced by Hastie and Tibshirani [13], which are defined as follows,

$$y=\sum _{j=1}^{p_n}{m}_j\left(x_j\right)+ϵ,$$

(1)

where y is the response variable, $x_j$ is the covariate, $m_j$ is an unknown function with $j=1,\dots ,p_n$, and $ϵ$ is the random error. Obviously, this additive combination of univariate functions could detect the nonlinear independence easily, but their good statistical properties and high computational performance only belong to low dimensions. For ultrahigh dimensions, to keep them working well, one of the most popular methods is the two-stage approach. Its main idea is to perform model selection in a fast and efficient way while retaining all the important features in the reduced feature space and then refitting the reduced models. In the following paper, we focus on the methodology of model selection for ultrahigh nonparametric additive models. In that field, the last decade has seen a growing trend toward smooth-group penalized methods, see [14,15,16,17]. Whereas the above methods may involve some tuning parameters, which bring a heavy computational burden and unstable results in high dimensions. A forward feature selection procedure, proposed by [11] for ultrahigh dimensional nonparametric additive models, does not involve any initial parameters. In addition, model-free methods have been developed recently. Based on the cumulative divergence (CD), Zhou et al. [12] proposed a forward screening procedure that considered the joint effects among covariates in the feature screening process. These two above methods screen the remaining candidate indexes into the sub-models through forwarding procedures. This kind of forward-searching algorithm also leads to a high computational burden. Furthermore, under previous studies’ correlation assumption, they ignored that the predictors are often correlated for high-dimensional feature space. In detail, the unimportant covariate $x_{\ell }$ corresponding to $m_{\ell }\equiv 0$ in the nonparametric additive models (1) may have a strong correlation with the residual $y-{\sum }_{j\in \mathcal{M}}{m}_j\left(x_j\right)$ given index set $\mathcal{M}\subseteq \{1,\dots ,p_n\}$, which implies that their methodologies may screen quite a few unimportant features into the sub-models.

To improve these limitations, first, our approach is to propose a group screening procedure via nonparametric smoothing ridge estimation (GRIE), motivated by the theoretical property and outstanding simulation performance of the ridge estimator in [18]. The core idea of GRIE is to get the importance of each covariate by combining the ridge estimator and group contribution. Its details are as follows. We begin with fitting the ridge regression by B-spline smoothing and then treating the spline basis corresponding to each covariate as a group. Next, we evaluate the group contribution of covariates by the magnitude of group estimators. Lastly, we sort the importance of the covariates by the group contribution in descending order. To further conduct model selection, we propose the refined GRIE-EBIC method mixing GRIE and the extended Bayesian information criteria (EBIC) in [19].  The GRIE-EBIC method is used to search for the predictor with the most group contributions by  EBIC.

Compared with other feature selection methods for nonparametric additive models, the GRIE-EBIC method has the following advantages: (1) the joint correlation among covariates is considered, and the strong marginal correlation assumption between response and important predictors is relaxed; (2) simple calculation with lower computational complexity; (3) strong consistency for feature screening, and it implies that the true features can be extracted accurately with probability tending to one, which does not exist in other stepwise feature screening methods, such as forward additive regression in [11], forward screening in [12], etc.

The rest of the paper is organized as follows. In Section 2, we introduce the GRIE screening procedure, the GRIE-EBIC method, and its algorithm. In Section 3, we establish the sure screening property of the GRIE screening procedure and the strong consistency of screening by the GRIE-EBIC. In Section 4, we present the performance of our proposed algorithm through simulation studies. In Section 5, we apply our methodology to fit two real datasets to further illustrate the performance of our proposed method. The first is based on Boston housing, while the second is related to Arabidopsis thaliana gene data. A conclusion is given in Section 6. The proofs are in Appendix A.

Notation 

Let $\mathbf{A}$ be m by l matrix, $\mathbf{M}$ be any subset of $\left\{1,2,\dots ,l\right\}$ with any positive integers of m and l, and then ${\mathbf{A}}_{\mathbf{M}}$ be the submatrix of $\mathbf{A}$ formed by column indexes in $\mathbf{M}$. We write ${\lambda }_{min}\left(\mathbf{A}\right)$ and ${\lambda }_{max}\left(\mathbf{A}\right)$ to denote the minimum and maximum eigenvalues of a symmetric matrix $\mathbf{A}$, separately. We write ${\mathbf{I}}_m$ as the identity matrix. We defined ${\mathbf{P}}_{\lambda ,\mathbf{A}}=\mathbf{A}{({\mathbf{A}}^{\top }\mathbf{A}+\lambda {\mathbf{I}}_m)}^{-1}{\mathbf{A}}^{\top }$, where λ is some positive constant, ${\mathbf{A}}^{\top }$ represents the transpose of matrix $\mathbf{A}$, and  here $\mathbf{A}$ is the column full rank $l\times m$ matrix with $m\le l$. When $\lambda =0$, ${\mathbf{P}}_{\mathbf{A}}=\mathbf{A}{\left({\mathbf{A}}^{\top }\mathbf{A}\right)}^{-1}{\mathbf{A}}^{\top }$, which is the projection onto column space of $\mathbf{A}$. Otherwise, ${\mathbf{e}}_i={(0,\dots ,0,1,0,\dots ,0)}^{\top }$ is the unit vector, which has zeros everywhere, except in the ith position. For vector $\mathbf{a}\in {\mathbb{R}}^n$, the $L^2$ norm of $\mathbf{a}={(a_1,a_2,\dots ,a_n)}^{\top }$ is captured by ${\parallel \mathbf{a}\parallel }_2=\sqrt{{\mathbf{a}}^{\top }\mathbf{a}}$.

2. Methodology

Suppose we have the random sample $\left\{(y_i,x_{i,1},\dots ,x_{i,p_n}):i=1,\dots ,n\right\}$, which is generated from the population model (1). Then the nonparametric additive model can be rewritten as:

$$y_i=\sum _{j=1}^{p_n}{m}_j\left(x_{i,j}\right)+{ϵ}_i,i=1,\dots ,n.$$

(2)

Without loss generality, we assume that the mean response is zero. For identification of the model, we further assume the mean of each additive function is zero, i.e., $E{m}_j\left(x_{i,j}\right)=0$ for $j=1,\dots ,p_n$. We note that all of the response variables are centralized to satisfy the above assumption during a real application. Here, the variance of the additive function $\mathrm{Var}\left(m_j\left(x_j\right)\right)$ is used to distinguish the importance of the covariate. Thus, we let $x_j$ be the important predictor if $\mathrm{Var}\left(m_j\left(x_j\right)\right)>0$; otherwise $x_j$ is the redundant predictor.  Then we define the index set of the important predictors as $\mathbf{S}=\{j:\mathrm{Var}\left(m_j\left(x_j\right)\right)>0,j=1,\dots ,p_n\}$.

Next, we use B-spline basis functions to approximate $m_j(·)$. Let us assume $x_j\in [0,1]$ for $j=1,\dots ,p_n$, $\overline{\varphi }={\left\{{\varphi }_k\right\}}_{k=0}^q$ be a knot sequence such that $0={\varphi }_0<{\varphi }_1<\dots <{\varphi }_q=1$, and  $\mathcal{S}(\ell ,\overline{\varphi })$ are the space of polynomial splines of order ℓ with knot sequence $\overline{\varphi }$. $\mathcal{S}(\ell ,\overline{\varphi })$ is a ${\kappa }_n$–dimensional linear space with ${\kappa }_n=q+\ell$. For any $m_j\left(x_j\right)$, $j=1,\dots ,p_n$, there exists the unique vector ${\theta }_j^{*}$ to satisfy

$$m_j\left(x_j\right)\approx \sum _{t=1}^{{\kappa }_n}{\theta }_{jt}{B}_t\left(x_j\right)=\mathbf{B}{\left(x_j\right)}^{\top }{\theta }_j^{*},$$

(3)

where $\mathbf{B}\left(x_j\right)={(B_1\left(x_j\right),\dots ,B_{{\kappa }_n}\left(x_j\right))}^{\top }$ and ${\theta }_j^{*}={({\theta }_{j1}^{*},\dots ,{\theta }_{j{\kappa }_n}^{*})}^{\top }$. Let ${\mathbf{w}}_i={({\mathbf{w}}_{i,1}^{\top },\dots ,{\mathbf{w}}_{i,p_n}^{\top })}^{\top }$ with ${\mathbf{w}}_{i,j}=\mathbf{B}\left(x_{i,j}\right)$, $\mathbf{W}={({\mathbf{w}}_1,\dots ,{\mathbf{w}}_n)}^{\top }$ and $\mathbf{Y}={(y_1,\dots ,y_n)}^{\top }$. Based on the approximation of (3), model (2) becomes

$$y_i={\mathbf{w}}_i^{\top }{\theta }^{*}+{ϵ}_i^{*},i=1,\dots ,n,$$

(4)

where ${\theta }^{*}={({{\theta }_1^{*}}^{\top },\cdots ,{{\theta }_{p_n}^{*}}^{\top })}^{\top }$ and ${ϵ}_i^{*}={\sum }_{j=1}^{p_n}{m}_j\left(x_{i,j}\right)-{\mathbf{w}}_i^{\top }{\theta }^{*}+{ϵ}_i$. Under model (4), the ridge estimator minimizes the following loss

$${\parallel \mathbf{Y}-\mathbf{W}\theta \parallel }_2^2+\lambda {\parallel \theta \parallel }_2^2,$$

where $\lambda$ is a positive constant. Then $\widehat{\theta }$ admits

$$\begin{array}{c}\hfill \widehat{\theta }={\mathbf{W}}^{\top }{(\mathbf{W}{\mathbf{W}}^{\top }+\lambda {\mathbf{I}}_n)}^{-1}\mathbf{Y},\end{array}$$

(5)

where ${\mathbf{I}}_n$ is the $n\times n$ identity matrix. For linear regression, Wang and Leng [18] considered the effect of each entry of $\theta$, and showed that the ridge estimator achieves screening consistency. Notice that $\mathrm{Var}\left(m_j\left(x_j\right)\right)\approx {{\theta }_j^{*}}^{\top }E\left({\mathbf{w}}_j{\mathbf{w}}_j^{\top }\right){\theta }_j^{*}$. Different from linear regression, we need to consider the group contribution of ${\theta }_j^{*}$. By the boundedness of $E\left({\mathbf{w}}_{i,j}{\mathbf{w}}_{i,j}^{\top }\right)$ from assumption A4(i), we use $\parallel {\theta }_j^{*}{\parallel }_2$ to evaluate the group contribution. Similar to the results in [18], the ridge estimator $\widehat{\theta }$ provides the ranking order of the group contribution in ${\theta }^{*}$ with $P\left(\parallel {\widehat{\theta }}_j{\parallel }_2>{\parallel {\widehat{\theta }}_k\parallel }_2\right)\to 1$ if $j\in \mathbf{S},k\in {\mathbf{S}}^c$ (see Theorem 1).

One natural screening method is to sort $\{\parallel {\widehat{\theta }}_j{\parallel }_2^2\}$ in decreasing order, and  select its top m indexes, denoted as ${\mathcal{F}}_m=\{i_1,i_2,i_3,\dots ,i_m\},1\le m\le p_n.$ This screening process is referred to as the “GRIE” screening procedure. We define $\mathcal{G}=\left\{{\mathcal{F}}_m:m=1,\dots ,p_n\right\}$, $\mathcal{A}=\left\{m:\mathbf{S}\subseteq {\mathcal{F}}_m,1\le m\le p_n\right\}.$ Further, to get a more accurate result of model selection, i.e., searching $d_n$, which is the minimum item in the set $\mathcal{A}$. At that time, ${\mathcal{F}}_{d_n}$ is a set with the shortest length in $\mathcal{G}$ that contains important variable set $\mathbf{S}$. With the definition of $\mathcal{G}$, we have $\mathbf{S}\subseteq {\mathcal{F}}_{p_n}$.  Then ${\mathcal{F}}_{d_n}$ is not an empty set. In summary, we want to find ${\mathcal{F}}_{d_n}$ from $\mathcal{G}$.

It is well known that the extended Bayesian information criteria (EBIC) have appealing theoretical properties and outstanding numerical performance for model selection. Let ${\mathbf{W}}_{\mathcal{T}}=({\mathbf{W}}_j,j\in \mathcal{T})$ for any subset of $\mathcal{T}\subset \{1,\dots ,p_n\}$. The formula of the EBIC for the sub-model $(\mathbf{Y},{\mathbf{W}}_{\mathcal{T}})$ is given by

$$EBIC\left(\mathcal{T}\right)=log(RSS\left(\mathcal{T}\right)/n)+\left\{{\kappa }_n\left|\mathcal{T}\right|log\left(n\right)+2\gamma logf\left(\right|\mathcal{T}\left|\right)\right\}/n,$$

(6)

where $\gamma$ is the preset positive constant, $RSS\left(\mathcal{T}\right)=\parallel \mathbf{Y}-{\mathbf{W}}_{\mathcal{T}}{\widehat{\theta }}_{\mathcal{T}}{\parallel }_2^2$ is the sum of squared residuals (RSS), and $f\left(\right|\mathcal{T}\left|\right)={\mathrm{C}}_{p_n{\kappa }_n}^{\left|\mathcal{T}\right|{\kappa }_n}$ is the combination number.

For a linear model, Wang [9] showed that $EBIC\left({\mathcal{F}}_m\right)<EBIC\left({\mathcal{F}}_{m-1}\right)$ if $i_m\in \mathbf{S}$. Based on this property of EBIC and the preserving rank property of GRIE screening procedure (see Theorem 1), we propose the following Algorithm 1 for the model selection of (1).

Algorithm 1: GRIE-EBIC algorithm.

Initialization: Input $(\mathbf{W}, \mathbf{Y})$, $RS{S}_0={\parallel \mathbf{Y}\parallel }_2^2$, n, $p_n$, $\lambda$, ${\kappa }_n$, $\gamma$, L. Step (i): Compute the GRIE screening procedure    1: Calculate ridge estimator $\widehat{\theta }={\mathbf{W}}^{\top }{({\mathbf{W}}^{\top }+\lambda {\mathbf{I}}_n)}^{-1}\mathbf{Y}$;    2: Sort $\{\parallel {\widehat{\theta }}_j{\parallel }_2,j=1,\dots ,p_n\}$ in decreasing order and select the top n index set    which is denoted by ${\mathcal{F}}_n=\{i_1,i_2,i_3,\dots ,i_n\}$; Step (ii): Direct decreasing solution path    3: For $k=1,\dots ,n$, do       $3.1$: Let ${\widehat{\mathbf{S}}}_k=\{i_1,\dots ,i_k\}$ and compute the sum of squared residuals                         $RS{S}_k={\parallel \mathbf{Y}-{\mathbf{W}}_{{\widehat{\mathbf{S}}}_k}{\left({\mathbf{W}}_{{\widehat{\mathbf{S}}}_k}^{\top }{\mathbf{W}}_{{\widehat{\mathbf{S}}}_k}\right)}^{-1}{\mathbf{W}}_{{\widehat{\mathbf{S}}}_k}^{\top }\mathbf{Y}\parallel }_2^2$ ;       $3.2$: Compute EBIC: $EBI{C}_k=log(RS{S}_k/n)+\{{\kappa }_{n}klog\left(n\right)+2\gamma logf\left(k\right)\}/n$;       $3.3$: If $k\ge L+1$ and $EBI{C}_k>\cdots >EBI{C}_{k-L}$, compute $K=k-L$ and stop;    4: Compute the difference of the EBIC to obtain the decreasing solution path                         $\mathcal{I}=\{k:EBI{C}_k-EBI{C}_{k-1}<0,k=1,2,\dots ,K\}$;    5: Find the decreasing index set ${\widehat{\mathbf{S}}}_{*}=\{i_k:k\in \mathcal{I}\}$; Step (iii): Forward decreasing solution path    6: Compute $RS{S}_{*}={\parallel \mathbf{Y}-{\mathbf{W}}_{{\widehat{\mathbf{S}}}_{*}}{\left({\mathbf{W}}_{{\widehat{\mathbf{S}}}_{*}}^{\top }{\mathbf{W}}_{{\widehat{\mathbf{S}}}_{*}}\right)}^{-1}{\mathbf{W}}_{{\widehat{\mathbf{S}}}_{*}}^{\top }\mathbf{Y}\parallel }_2^2$ and                         $EBI{C}_{*}=log(RS{S}_{*}/n)+\{{\kappa }_n|{\widehat{\mathbf{S}}}_{*}|log\left(n\right)+2\gamma logf\left(\right|{\widehat{\mathbf{S}}}_{*}\left|\right)\}/n;$    7: For $\ell \in {\mathcal{F}}_n\setminus {\widehat{\mathbf{S}}}_{*}$, do       Let ${\widehat{\mathbf{S}}}_{*\ell }={\widehat{\mathbf{S}}}_{*}\cup \{\ell \}$, compute $RS{S}_{*\ell }={\parallel \mathbf{Y}-{\mathbf{W}}_{{\widehat{\mathbf{S}}}_{*\ell }}{\left({\mathbf{W}}_{{\widehat{\mathbf{S}}}_{*\ell }}^{\top }{\mathbf{W}}_{{\widehat{\mathbf{S}}}_{*\ell }}\right)}^{-1}{\mathbf{W}}_{{\widehat{\mathbf{S}}}_{*\ell }}^{\top }\mathbf{Y}\parallel }_2^2$ and                         $EBI{C}_{*\ell }=log(RS{S}_{*\ell }/n)+\{{\kappa }_n|{\widehat{\mathbf{S}}}_{*\ell }|log\left(n\right)+2\gamma logf\left(\right|{\widehat{\mathbf{S}}}_{*\ell }\left|\right)\}/n;$    8: Find decreasing solution path $\widehat{\mathbf{S}}={\widehat{\mathbf{S}}}_{*}\cup \{\ell :EBI{C}_{*\ell }-EBI{C}_{*}<0,\ell \in {\mathcal{F}}_n\setminus {\widehat{\mathbf{S}}}_{*}\}$; Output final index set $\widehat{\mathbf{S}}$.

In step (ii) of the GRIE-EBIC algorithm, we search the important covariates from the top n predictor space ${\mathcal{F}}_n$. Based on Theorem 1, GRIE has the consistency of preserving order in sorting. The higher the index position of the variable in ${\mathcal{F}}_{p_n}$, the more likely it is to be an important variable. To speed up the calculation, we set a stopping rule for screening when the EBIC value increases for L times continuously. To improve the robustness of the GRIE-EBIC algorithm, in step (iii), we add the further forward screening process.

3. Asymptotic Properties

3.1. Assumptions

To establish the asymptotic properties of our proposed method, we give the following notations and assumptions. Let $\mathsf{\Sigma }=E\left(\mathbf{w}{\mathbf{w}}^{\top }\right)$, $\mathbf{Z}=\mathbf{W}{\mathsf{\Sigma }}^{-1/2}$, $\mathbf{z}={\mathsf{\Sigma }}^{-1/2}\mathbf{w}$, and $t_n=p_n{\kappa }_n$, where $\mathbf{w}={({\mathbf{w}}_1^{\top },\dots ,{\mathbf{w}}_{p_n}^{\top })}^{\top }$ with ${\mathbf{w}}_j=\mathbf{B}\left(x_j\right)$. We use ${\mathsf{\Sigma }}_{\mathcal{T}}=E\left(n^{-1}{\mathbf{W}}_{\mathcal{T}}^{\top }{\mathbf{W}}_{\mathcal{T}}\right)$. ${\mathcal{H}}_r$ denotes a space of functions whose d-th order derivative is H$\ddot{o}$lder continuous of order v, i.e., ${\mathcal{H}}_r=\{h\left(z\right):|h^{\left(d\right)}\left(a^{\prime }\right)-h^{\left(d\right)}\left(a\right)|\le C|a^{\prime }-a{|}^v,a,a^{\prime }\in [0,1]\}$, where $h^{\left(d\right)}(·)$ is the d-th derivative of $h(·)$ and $r=d+v$. If $v=1$, $h^{\left(d\right)}(·)$ is Lipschitz continuous. Let $s_n$ be the cardinality of $\mathbf{S}$. The following assumptions are required:

A1

Assume $\mathbf{z}$ has a spherically symmetric distribution, and there exists some positive $c_1$ and $C_1$ such that

$$P\left({\lambda }_{min}\left({t_n}^{-1}\mathbf{Z}{\mathbf{Z}}^{\top }\right)\le c_1^{-1}\mathrm{or}{\lambda }_{max}\left({t_n}^{-1}\mathbf{Z}{\mathbf{Z}}^{\top }\right)>c_1\right)\le 2exp(-C_{1}n).$$

A2

Assume there exists some positive constant $C_{*}$ such that, for any $a\in R$,

$$\underset{i=1,\dots ,n}{max}E\left\{exp\left(a{\epsilon }_i\right)|{\mathbf{x}}_i\right\}\le exp(C_{*}{a}^2/2).$$

A3

Assume that (i) there exists some $r\ge 2$ such that $m_j\in {\mathcal{H}}_r$ and ${\kappa }_n=O\left(n^{1/(2r+1)}\right)$ for any $j\in \mathbf{S}$; (ii) ${\sum }_{j\in \mathbf{S}}E{\left|m_j\left(x_j\right)\right|}^2\le c_2{s}_n$; (iii) ${\lambda }_{max}\left(\mathsf{\Sigma }\right)/{\lambda }_{min}\left(\mathsf{\Sigma }\right)\le c_3{n}^{\tau }$, where $c_2,c_3$ are some positive constants and $\tau \ge 0$.

A4

(i) $c_4^{-1}{\kappa }_n^{-1}\le {\lambda }_{min}\left(E\left(\mathbf{B}\left(x_j\right)\mathbf{B}{\left(x_j\right)}^{\top }\right)\right)\le {\lambda }_{max}\left(E\left(\mathbf{B}\left(x_j\right)\mathbf{B}{\left(x_j\right)}^{\top }\right)\right)\le c_4{\kappa }_n^{-1}$ for some positive constant $c_4$; (ii) ${min}_{j\in \mathbf{S}}\left\{E\right|m_j\left(x_j\right){{|}^2\}}^{1/2}\ge d_n$ for some positive sequence $d_n\to 0$; (iii) $\frac{{\kappa }_n^{r-1/2}{d}_n}{n^{2\tau }{s}_n\sqrt{logn}}\to \infty$, $log\left(t_n\right)=o\left(\frac{d_n^2{n}^{1-4\tau }}{{\kappa }_n^2{s}_n^{2}logn}\right)$.

A5

(i) $\mathrm{Var}\left(y_1\right)=O({\kappa }_n{s}_n^2{n}^{3\tau }log\left(n\right))$; (ii) For any integer N with $s_n<N\le s_{n}logn$, there exists positive constant $c_6>0$ such that

$$\frac{c_6{n}^{-\tau }}{{\kappa }_n}\le {\lambda }_{min}\left({\mathsf{\Sigma }}_{\mathcal{T}}\right)$$

(7)

holds uniformly in $\mathcal{T}\subset {\mathcal{F}}_n$ satisfying $\left|\mathcal{T}\right|\le N$ and $\mathbf{S}\subset \mathcal{T}$.

Assumptions A1 and A3(iii) are like Assumptions 1 and 3 of [18]. Assumption A2 is the same as Assumption A3 of [11], which means that the random error follows the sub-Gaussian distribution. Assumption A3(i) is a common assumption in the literature for the polynomial spline basis, A3(ii) gives the upper bound of all signals, and A3(iii) gives the upper bound of the condition number. In addition, Assumption A3(ii)–(iii) are implied by Assumption A2 in [11]. Assumption A4(i) and a stronger assumption $\mathrm{Var}\left(y_1\right)=O\left(1\right)$ than A5(i), which is also imposed in [11] for achieving the consistency of variable selection. They also assumed A5(ii) holds. Assumptions A4(ii) and (iii) give the lower bound and upper bound of the minimal signal and dimensionality of the design matrix $\mathbf{W}$.

3.2. Main Theorems

Theorem 1.

If Assumptions A1–A4 hold, then

$$P\left(\underset{j\in \mathbf{S}}{min}\parallel {\widehat{\theta }}_j{\parallel }_2>\underset{j\in {\mathbf{S}}^c}{max}{\parallel {\widehat{\theta }}_j\parallel }_2\right)\to 1.$$

Alternatively, we can choose a sub-model ${\mathcal{F}}_{d_n}$ with $d_n=O\left(n^{\iota }\right)$ for some $0<\iota <1$ such that

$$P\left(\mathbf{S}\subset {\mathcal{F}}_{d_n}\right)\to 1.$$

Theorem 1 states the consistency of preserving order in sorting, i.e., $\widehat{\theta }$ could totally separate the unimportant and important variables with a probability tending to 1. For the linear models, Theorem 1 is in line with Theorem 2 in [18], which is the special case of our theorem.

Theorem 2.

If Assumptions A1–A5 hold, then

$$P\left(\widehat{\mathbf{S}}=\mathbf{S}\right)\to 1.$$

The screening methods in [7,11,12] adopted a forward selection algorithm, which means the later results are affected by the results of the previous steps. This not only brings a heavy computation burden but also results in overfitting results for screening with $P(\mathbf{S}\subseteq \widehat{\mathbf{S}})\to 1$. Compared with this result, Theorem 2 gives strong consistency of screening with $P\left(\widehat{\mathbf{S}}=\mathbf{S}\right)\to 1$.

4. Simulations

In this section, we investigate the finite sample performance of our proposed method and compare our method with the following two procedures: forward additive regression (FAR) in [11] and cumulative divergence-based forward regression (C-FS) in [12]. We choose $\lambda =1$, $L=5$ (suggested by [10]), $\gamma =0.5$ (suggested by [20]), and ${\kappa }_n=\lfloor n^{1/5}\rfloor +2$ (suggested by [11]) for the GRIE-EBIC algorithm, where ${\kappa }_n$ is the dimension of B-spline basis space and $\lfloor n^{1/5}\rfloor$ is the greatest integer less than $n^{1/5}$.

Three specific criteria are adopted to evaluate the performance of variable selection for the additive model (1). True positive (TP) is the number of the true variables that are considered true variables in the selected model, and false positive (FP) is the number of the noise variables that are misclassified as true variables in the selected model. Combining TP and FP, they reflect the accuracy of variable selection methods in the selected sub-models. In addition, we select time as the third criterion to reflect the efficiency of variable selection by different methods. It is easy to find that our proposed method, GRIE, is more efficient than FAR and C-FS in computation since their complexities of calculation are $O\left(n^2{p}_n{\kappa }_n\right)$, $O\left(n^3{p}_n\right)$, and $O\left(T{n}^3{p}_n\right)$, respectively, where T is the number of repetitions for the bootstrap procedure in the C-FS method. The comparison of computation complexities highlights the time efficiency of GRIE in the calculation, which will be further demonstrated by simulation results in

Table 1. Average number of true positive (TP), false positive (FP), and calculation time over 100 repetitions and their robust standard deviations (in parentheses) of Example 1 with $ϵ\sim N(0,1)$.

$\mathit{\rho }$	Approach	${\mathit{p}}_{\mathit{n}}=500$			${\mathit{p}}_{\mathit{n}}=1000$
		TP	FP	Time (s)	TP	FP	Time (s)
AR Structure
0.3	FAR	4.00 (0.00)	0.00 (0.00)	83.19 (9.80)	4.00 (0.00)	0.00 (0.00)	166.26 (18.65)
	C-FS	3.20 (0.40)	5.21 (2.96)	16.18 (5.38)	3.34 (0.48)	11.44 (5.38)	39.64 (14.35)
	GRIE	4.00 (0.00)	0.00 (0.00)	2.37 (0.28)	3.99 (0.10)	0.01 (0.10)	3.56 (0.77)
0.6	FAR	4.00 (0.00)	0.00 (0.00)	82.06 (9.77)	4.00 (0.00)	0.00 (0.00)	168.16 (20.51)
	C-FS	3.71 (0.46)	4.81 (2.39)	16.57 (4.28)	3.70 (0.46)	9.33 (4.88)	34.61 (12.24)
	GRIE	3.99 (0.10)	0.00 (0.00)	2.40 (0.35)	3.98 (0.14)	0.03 (0.30)	3.43 (0.72)
0.9	FAR	3.17 (0.60)	0.00 (0.00)	81.34 (9.48)	3.09 (0.60)	0.00 (0.00)	168.70 (18.62)
	C-FS	3.14 (0.51)	2.44 (1.72)	10.63 (3.01)	3.14 (0.53)	4.43 (2.96)	19.14 (6.78)
	GRIE	3.71 (0.46)	0.20 (0.40)	2.22 (0.40)	3.70 (0.46)	0.21 (0.43)	3.45 (0.76)
CS Structure
0.3	FAR	4.00 (0.00)	0.00 (0.00)	83.60 (10.17)	4.00 (0.00)	0.00 (0.00)	165.38 (19.00)
	C-FS	3.45 (0.52)	4.96 (2.97)	16.11 (5.23)	3.33 (0.47)	11.69 (6.24)	39.98 (16.38)
	GRIE	4.00 (0.00)	0.09 (0.90)	2.30 (0.39)	4.00 (0.00)	0.02 (0.20)	3.57 (0.72)
0.6	FAR	4.00 (0.00)	0.00 (0.00)	84.24 (10.26)	4.00 (0.00)	0.01 (0.10)	166.92 (18.64)
	C-FS	3.74 (0.44)	5.05 (2.98)	16.82 (5.26)	3.61 (0.55)	10.26 (5.25)	36.73 (13.38)
	GRIE	4.00 (0.00)	0.23 (2.30)	2.35 (0.37)	4.00 (0.00)	0.11 (0.62)	3.41 (0.74)
0.9	FAR	3.03 (0.67)	0.00 (0.00)	85.57 (11.01)	2.79 (0.70)	0.00 (0.00)	166.02 (18.65)
	C-FS	2.63 (0.65)	4.48 (3.47)	13.35 (6.08)	2.56 (0.67)	9.70 (5.93)	32.23 (15.07)
	GRIE	3.89 (0.31)	1.47 (7.62)	2.21 (0.35)	3.79 (0.41)	3.16 (17.90)	3.44 (0.77)

and

Table 2. Average numbers of true positive (TP), false positive (FP), and calculation time over 100 repetitions and their robust standard deviations (in parentheses) of Example 1 with $ϵ\sim 0.5{\chi }_2^2$.

$\mathit{\rho }$	Approach	${\mathit{p}}_{\mathit{n}}=500$			${\mathit{p}}_{\mathit{n}}=1000$
		TP	FP	Time (s)	TP	FP	Time (s)
AR Structure
0.3	FAR	4.00 (0.00)	0.00 (0.00)	79.30 (11.14)	4.00 (0.00)	0.00 (0.00)	165.38 (22.09)
	C-FS	3.27 (0.45)	5.38 (3.00)	16.43 (5.40)	3.33 (0.47)	11.24 (4.85)	39.44 (12.59)
	GRIE	4.00 (0.00)	0.00 (0.00)	2.40 (0.36)	3.99 (0.10)	0.00 (0.00)	3.33 (0.68)
0.6	FAR	4.00 (0.00)	0.00 (0.00)	79.19 (12.07)	4.00 (0.00)	0.00 (0.00)	163.64 (23.37)
	C-FS	3.70 (0.46)	4.42 (2.53)	15.60 (4.49)	3.77 (0.42)	9.56 (4.29)	35.74 (11.14)
	GRIE	3.99 (0.10)	0.01 (0.10)	2.31 (0.33)	3.98 (0.14)	0.03 (0.30)	3.42 (0.67)
0.9	FAR	3.09 (0.68)	0.00 (0.00)	80.28 (10.89)	3.01 (0.72)	0.00 (0.00)	163.88 (22.78)
	C-FS	3.10 (0.48)	2.28 (1.56)	10.15 (2.86)	3.15 (0.50)	4.26 (2.20)	19.12 (5.44)
	GRIE	3.71 (0.46)	0.23 (0.51)	2.28 (0.32)	3.78 (0.42)	0.16 (0.39)	3.33 (0.67)
CS Structure
0.3	FAR	4.00 (0.00)	0.00 (0.00)	80.28 (9.59)	4.00 (0.00)	0.00 (0.00)	164.25 (19.31)
	C-FS	3.51 (0.52)	5.10 (2.99)	16.50 (5.68)	3.36 (0.48)	10.87 (4.98)	37.91 (13.02)
	GRIE	4.00 (0.00)	0.00 (0.00)	2.31 (0.34)	3.98 (0.14)	0.34 (3.30)	3.38 (0.66)
0.6	FAR	4.00 (0.00)	0.00 (0.00)	80.12 (11.56)	4.00 (0.00)	0.01 (0.10)	165.41 (20.04)
	C-FS	3.72 (0.49)	4.71 (2.57)	16.04 (4.68)	3.68 (0.49)	9.79 (5.39)	36.12 (14.41)
	GRIE	3.99 (0.10)	0.00 (0.00)	2.31 (0.29)	4.00 (0.00)	0.11 (0.65)	3.40 (0.70)
0.9	FAR	3.00 (0.79)	0.03 (0.17)	79.62 (11.60)	2.85 (0.78)	0.02 (0.14)	164.90 (19.25)
	C-FS	2.73 (0.66)	4.40 (2.81)	13.56 (4.91)	2.69 (0.63)	10.82 (5.99)	35.72 (15.68)
	GRIE	3.94 (0.24)	3.13 (16.94)	2.28 (0.35)	3.82 (0.39)	3.63 (18.13)	3.30 (0.72)

.

The following examples perform the effect from different dimensions and correlations between any two covariates with each other by the above three procedures. Given two different dimensions and three different correlations between any two predictors with each other, the considered error followed the standard normal $N(0,1)$ and Chi-square $0.5{\chi }_2^2$. In each example, we generate 100 random samples, each consisting of $n=300$. The data generation procedures are implemented by the R package “MASS” generating the random covariates, errors, and response variables. The details are as follows: (1) “rnorm”: simulated from a multivariate normal distribution; (2) “mvrnorm”: simulated from a multivariate normal distribution; (3) “rchisqure”: simulated from a multi Chi-square distribution; (4) “runif”: simulated from a uniform distribution.

Example 1.

We generated n samples from the following nonparametric additive model:

$$y=m_1\left(x_1\right)+m_2\left(x_2\right)+m_3\left(x_3\right)+m_4\left(x_4\right)+ϵ,$$

where $m_1\left(x\right)=0.75exp\left(x\right)$, $m_2\left(x\right)=x^2$, $m_3\left(x\right)=3sin\left(x\right)$, $m_4\left(x\right)=2x$, and ${(x_1,x_2,\dots ,x_{p_n})}^{\top }$ follows a multinormal distribution $N(0,\mathsf{\Sigma })$. In this example, given $\mathsf{\Sigma }=\left({\sigma }_{ij}\right)$ under the following two cases: $\left(1\right)$ Autoregressive (AR) structure, ${\sigma }_{ij}={\rho }^{|i-j|}$; $\left(2\right)$ Compound symmetry (CS) structure, namely, if $i\ne j$, ${\sigma }_{ij}=\rho$, else ${\sigma }_{ij}=1$. Here, we set the parameter ρ used to control the strength of correlation between any two predictors with each other at different values of $0.3$, $0.6$, and $0.9$.

Table 1, Table 2,

Table 3. The empirical probabilities of each important covariate and all important covariates being retained for 100 replications in Example 1 with $ϵ\sim N(0,1)$.

$\mathit{\rho }$	Approach	${\mathit{p}}_{\mathit{n}}=500$					${\mathit{p}}_{\mathit{n}}=1000$
		${\mathcal{P}}_{\mathbf{1}}$	${\mathcal{P}}_{\mathbf{2}}$	${\mathcal{P}}_{\mathbf{3}}$	${\mathcal{P}}_{\mathbf{4}}$	${\mathcal{P}}_{\mathbf{all}}$	${\mathcal{P}}_{\mathbf{1}}$	${\mathcal{P}}_{\mathbf{2}}$	${\mathcal{P}}_{\mathbf{3}}$	${\mathcal{P}}_{\mathbf{4}}$	${\mathcal{P}}_{\mathbf{all}}$
AR Structure
0.3	FAR	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	C-FS	1.00	0.20	1.00	1.00	0.20	1.00	0.34	1.00	1.00	0.34
	GRIE	1.00	1.00	1.00	1.00	1.00	1.00	0.99	1.00	1.00	0.99
0.6	FAR	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	C-FS	1.00	0.71	1.00	1.00	0.71	0.98	0.72	1.00	1.00	0.70
	GRIE	1.00	0.99	1.00	1.00	0.99	1.00	0.98	1.00	1.00	0.98
0.9	FAR	0.80	0.49	0.97	0.91	0.28	0.82	0.42	0.98	0.87	0.23
	C-FS	0.50	0.67	0.97	1.00	0.21	0.50	0.67	0.97	1.00	0.22
	GRIE	0.83	0.88	1.00	1.00	0.71	0.79	0.91	1.00	1.00	0.70
CS Structure
0.3	FAR	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	C-FS	0.99	0.46	1.00	1.00	0.46	1.00	0.33	1.00	1.00	0.33
	GRIE	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
0.6	FAR	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	C-FS	0.93	0.81	1.00	1.00	0.74	0.88	0.73	1.00	1.00	0.64
	GRIE	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
0.9	FAR	0.63	0.57	0.98	0.85	0.24	0.58	0.49	0.91	0.81	0.16
	C-FS	0.11	0.62	0.97	0.93	0.05	0.12	0.53	0.97	0.94	0.07
	GRIE	0.93	0.96	1.00	1.00	0.89	0.87	0.92	1.00	1.00	0.79

and

Table 4. The empirical probabilities of each important covariate and all important covariates being retained for 100 replications in Example 1 with $ϵ\sim 0.5{\chi }_2^2$.

$\mathit{\rho }$	Approach	${\mathit{p}}_{\mathit{n}}=500$					${\mathit{p}}_{\mathit{n}}=1000$
		${\mathcal{P}}_{\mathbf{1}}$	${\mathcal{P}}_{\mathbf{2}}$	${\mathcal{P}}_{\mathbf{3}}$	${\mathcal{P}}_{\mathbf{4}}$	${\mathcal{P}}_{\mathbf{all}}$	${\mathcal{P}}_{\mathbf{1}}$	${\mathcal{P}}_{\mathbf{2}}$	${\mathcal{P}}_{\mathbf{3}}$	${\mathcal{P}}_{\mathbf{4}}$	${\mathcal{P}}_{\mathbf{all}}$
AR Structure
0.3	FAR	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	C-FS	1.00	0.27	1.00	1.00	0.27	1.00	0.33	1.00	1.00	0.33
	GRIE	1.00	1.00	1.00	1.00	1.00	1.00	0.99	1.00	1.00	0.99
0.6	FAR	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	C-FS	1.00	0.70	1.00	1.00	0.70	0.99	0.78	1.00	1.00	0.77
	GRIE	1.00	0.99	1.00	1.00	0.99	1.00	0.99	1.00	0.99	0.98
0.9	FAR	0.81	0.47	0.95	0.86	0.28	0.82	0.45	0.94	0.80	0.26
	C-FS	0.42	0.70	0.98	1.00	0.17	0.53	0.65	0.97	1.00	0.21
	GRIE	0.82	0.89	1.00	1.00	0.71	0.84	0.94	1.00	1.00	0.78
CS Structure
0.3	FAR	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	C-FS	0.99	0.52	1.00	1.00	0.52	0.99	0.37	1.00	1.00	0.36
	GRIE	1.00	1.00	1.00	1.00	1.00	1.00	0.98	1.00	1.00	0.98
0.6	FAR	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	C-FS	0.93	0.79	1.00	1.00	0.74	0.92	0.76	1.00	1.00	0.69
	GRIE	1.00	0.99	1.00	1.00	0.99	1.00	1.00	1.00	1.00	1.00
0.9	FAR	0.54	0.65	0.97	0.84	0.28	0.55	0.55	0.94	0.81	0.24
	C-FS	0.12	0.66	0.97	0.98	0.09	0.14	0.61	0.97	0.97	0.08
	GRIE	0.97	0.97	1.00	1.00	0.94	0.88	0.94	1.00	1.00	0.82

summarize the results for the additive model in Example 1. Under the setting of $\rho =0.3$ and $0.6$, except for C-FS, our proposed method and FAR method could identify all important features and keep the FP value close to zero in both AR and CS structures. Even so, the FAR method has the longest calculation time among the three methods. Furthermore, when there exists strong correlations between covariates ($\rho =0.9$), the performances of all three methods are worse at identifying important variables, especially for FAR and C-FS methods. Under this situation, compared with the other two methods, our method has the highest TP and shortest cost time. To perform the stability of our model, we report the empirical probabilities of each important covariate and all important covariates are retained for 100 replications in Table 3 and Table 4, where ${\mathcal{P}}_j$ and ${\mathcal{P}}_{\mathrm{all}}$ are the empirical probabilities of each important covariate and all important covariates being retained in the selected sub-model, respectively. Following Table 3 and Table 4, ${\mathcal{P}}_{\mathrm{all}}$ is below $0.3$ for FAR and C-FS, while the ${\mathcal{P}}_{\mathrm{all}}$ of GRIE is at least $0.70$. In addition, the ${\mathcal{P}}_{\mathrm{j}}$’s of our method is the best among the three methods in high-dimensional settings. Hence, we conclude that our proposed GRIE method performs robustly in the model selection of nonparametric additive models under high-dimension settings.

Example 2.

In this example, we consider a linear model with a group structure given by

$$y=\sum _{i=1}^{p_n}{\beta }_i{x}_i+ϵ$$

with the predictors being generated by the following process

$$\left\{\begin{array}{cc}{x}_i=z_1+z+w_i,\hfill & \forall i=1,3,\\ x_i=z_2+z+w_i,\hfill & \forall i=2,4,\\ x_5,\dots ,x_{p_n}\stackrel{i.i.d}{\sim }N(0,1),\hfill \end{array}\right.$$

where $w_1,\dots ,w_4\stackrel{i.i.d}{\sim }U(0,1)$, $z_1,z_2\stackrel{i.i.d}{\sim }$ $U(0,1)$ and the common component $z\sim N\left(0,{\delta }^2\right)$. The variance parameter δ is set at different values of $0.4$, $0.6$, and $0.8$ to control the strength of the group structure. The true value of the coefficients are ${\beta }_i=3$ with $i=1,\dots ,4$ and ${\beta }_i=0$ with $i=5,\dots ,p_n$.

We also conducted simulations with the normal errors and chi-square errors for Example 2 and found that the performances of the two errors in this example were very close. Therefore, we omitted the results of chi-square errors to save space. In the following, we only report the results from the normal errors in

Table 5. Average numbers of true positive (TP), false positive (FP), and calculation time over 100 repetitions and their robust standard deviations (in parentheses) of Example 2 with $ϵ\sim N(0,1)$.

$\mathit{\delta }$	Approach	${\mathit{p}}_{\mathit{n}}=500$			${\mathit{p}}_{\mathit{n}}=1000$
		TP	FP	Time (s)	TP	FP	Time (s)
0.4	FAR	4.00 (0.00)	0.59 (0.51)	81.85 (10.28)	3.98 (0.20)	0.57 (0.50)	168.16 (20.16)
	C-FS	4.00 (0.00)	5.32 (2.97)	18.74 (5.35)	4.00 (0.00)	11.40 (5.37)	42.87 (15.18)
	GRIE	4.00 (0.00)	0.04 (0.24)	2.41 (0.36)	4.00 (0.00)	0.06 (0.34)	3.59 (0.64)
0.6	FAR	3.94 (0.28)	1.09 (0.49)	80.25 (9.06)	3.86 (0.49)	1.07 (0.48)	164.35 (18.90)
	C-FS	4.00 (0.00)	6.05 (2.88)	19.32 (5.19)	4.00 (0.00)	12.11 (5.37)	43.61 (14.17)
	GRIE	4.00 (0.00)	0.17 (0.49)	2.33 (0.34)	4.00 (0.00)	0.18 (0.54)	3.43 (0.63)
0.8	FAR	3.66 (0.73)	1.26 (0.50)	80.04 (9.10)	3.68 (0.72)	1.22 (0.54)	164.05 (18.76)
	C-FS	3.81 (0.42)	5.88 (2.82)	18.62 (5.41)	3.85 (0.36)	12.13 (5.18)	42.92 (13.12)
	GRIE	3.95 (0.22)	0.38 (0.72)	2.37 (0.34)	3.89 (0.31)	0.27 (0.63)	3.45 (0.75)

and

Table 6. The empirical probabilities of each important covariate and all important covariates being retained for 100 replications in Example 2 with $ϵ\sim N(0,1)$.

$\mathit{\delta }$	Approach	${\mathit{p}}_{\mathit{n}}=500$					${\mathit{p}}_{\mathit{n}}=1000$
		${\mathcal{P}}_{\mathbf{1}}$	${\mathcal{P}}_{\mathbf{2}}$	${\mathcal{P}}_{\mathbf{3}}$	${\mathcal{P}}_{\mathbf{4}}$	${\mathcal{P}}_{\mathbf{all}}$	${\mathcal{P}}_{\mathbf{1}}$	${\mathcal{P}}_{\mathbf{2}}$	${\mathcal{P}}_{\mathbf{3}}$	${\mathcal{P}}_{\mathbf{4}}$	${\mathcal{P}}_{\mathbf{all}}$
0.4	FAR	1.00	1.00	1.00	1.00	1.00	1.00	1.00	0.99	0.99	0.99
	C-FS	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	GRIE	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
0.6	FAR	0.99	0.99	0.98	0.98	0.95	0.96	0.98	0.97	0.95	0.92
	C-FS	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	GRIE	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
0.8	FAR	0.92	0.91	0.91	0.92	0.81	0.90	0.92	0.93	0.93	0.83
	C-FS	0.93	0.97	0.96	0.95	0.82	0.94	0.97	0.97	0.97	0.85
	GRIE	1.00	1.00	0.97	0.98	0.95	0.97	1.00	0.95	0.97	0.89

. We find that the FAR’s performance is the worst for identifying important features with the increase in correlations between groups, the performances of C-FS and our method GRIE also become worse when $\delta$ is over $0.6$, while GRIE performs better even if there exists a strong correlation among covariates. The above phenomena can be further explained by Table 6. When $\delta$ is over $0.6$, there is no longer an overwhelming empirical probability of screening important covariates for FAR and C-FS, which results in a decrease in TP and ${\mathcal{P}}_{\mathrm{all}}$ values. However, our proposed method is still relatively robust to different values of $\delta$ in terms of TP and ${\mathcal{P}}_{\mathrm{all}}$.

5. Real Data

5.1. Boston Housing Data

We use the Boston housing dataset to further illustrate the performance of our proposed method. The dataset contains the MEDV (median value of owner-occupied homes) in 506 U.S. census tracts of Boston from the 1970 census and 13 other variables that explain the variation in housing value. The 13 explaining variables are RM (average number of rooms per dwelling), AGE (proportion of owner-occupied units built prior to 1940), RAD (index of accessibility to radial highways), TAX (full-value property-tax rate per 10,000), PTRATIO (pupil-teacher ratio by town), B ($1000{(\mathrm{Bk}-0.63)}^2$, Bk is the proportion of blacks by town), LSTAT (lower status of the population), CRIM (per capita crime rate by town), ZN (proportion of residential land zoned for lots over 25,000 square footage), INDUS (proportion of non-retail business acres per town), CHAS (Charles River dummy variable), NOX (nitric oxides concentration parts per 10 million), and DIS (weighted distances to five Boston employment centers). To simplify notation, we denote the covariates RM, AGE, RAD, TAX, PTRATIO, B, LSTAT, CRIM, ZN, INDUS, CHAS, NOX, and DIS as $x_1,\dots ,x_{13}$. To study the relationship between MEDV and the above 13 variables, we consider the following nonparametric addition models:

$$y=\sum _{j=1}^{13}{m}_j\left(x_j\right)+ϵ,$$

(8)

where y is the log(MEDV). In order to extend the above model to the setting of high-dimensional data, followed by [21], we generate artificial variables $x_j$ to add noise variables, which is defined as follows.

$$x_j=\frac{Z_j+2W}{3},$$

for $j=14,\dots ,1000$ into (8), where $Z_{14},\dots ,Z_{1000}\stackrel{i.i.d}{\sim }N(0,1)$, and $W\sim U(0,1)$.

We use FAR, C-FS, and GRIE to identify important variables in the above additive model (8) with the full dataset. The results are as follows.

(i)

Under FAR, 3 covariates $\{x_1,x_7,x_8\}$ are selected, denoted by “model ($A_1$)”.

(ii)

Under GRIE, we receive 6 covariates $\{x_1,x_5,x_6,x_7,x_8,x_{12}\}$, denoted by “model ($B_1$)”.

(iii)

Under C-FS, there are 15 covariates chosen. They are $\{x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8,$ $x_9,x_{12},x_{13},x_{156},x_{377},x_{737},x_{859}\}$, denoted by “model ($C_1$)”.

The above three sub-models have such a nest relation with $A_1\subseteq B_1\subseteq C_1$, and we want to investigate which model is best for fitting this dataset. The nondegenerative Vuong test in [22] is considered here to compare two nested models, and its null hypothesis is that the two models are equivalent. We first compare model ($A_1$) with model ($B_1$) by the Vuong test, and its p-value $=0.001$, which leads to the rejection of the null hypothesis. This indicates model ($B_1$) is better than model ($A_1$) since ($A_1$) is nested in ($B_1$). We also compare model ($B_1$) with model ($C_1$), and its corresponding p-value of the Vuong test equals $0.981$, which indicates models ($B_1$) and ($C_1$) are equivalent since the null hypothesis is not rejected in this situation. However, model ($B_1$) has a smaller model size than model ($C_1$). Therefore, model ($B_1$) is more suitable than model ($C_1$) for fitting the Boston housing dataset to be the best working model, which indicates that GRIE performs the best in identifying the important variables among the above three variable selection methods.

To further demonstrate our results, we compare FAR, C-FS, and GRIE through their prediction errors. Toward this end, we randomly select 100 validation sets, with each of which the full sample is randomly partitioned into the training and validation sets with the size ratio $4:1$. The training sets are for variable section, and the validation sets are for the estimation of the prediction error. We centralize the response variable y and choose the cubic splines ${\kappa }_n=3$ to approximate the additive function. The average numbers of model size, the number of selected noise variables (SNV), and adjusted mean prediction errors (A-PE) are used to evaluate the performance of the three methods. All the results are reported in

Table 7. Average numbers of model size, the number of SNV, and A-PE over 100 repetitions and their robust standard deviations (in parentheses) of Boston Housing Data.

Approach	Model Size	SNV	A-PE
FAR	2.10 (0.30)	0.00 (0.00)	0.052 (0.011)
C-FS	19.26 (5.39)	8.71 (5.10)	0.047 (0.012)
GRIE	5.07 (0.95)	0.00 (0.00)	0.043 (0.010)

.

From Table 7, we have that: $\left(1\right)$ The model sizes of our method, GRIE, and FAR are both smaller than C-FS, but the A-PE of FAR is the largest among the three methods, which means that FAR may fail to identify some important variables. To verify it, we report the frequency for 13 real covariates being selected over 100 replications in

Table 8. The frequency for 13 real covariates being selected over 100 replications for Boston Housing Data.

Variable	FAR	C-FS	GRIE
RM	100	100	100
AGE	0	99	0
RAD	0	60	6
TAX	0	59	7
PTRATIO	0	100	68
B	0	92	99
LSTAT	100	100	100
CRIM	10	100	80
ZN	0	97	0
INDUS	0	22	0
CHAS	0	26	0
NOX	0	100	47
DIS	0	100	0

. Table 8 shows that RM and LSTAT are selected by all methods in each repetition. Except for the FAR method, PTRATIO, B, and CRIM can be selected by GRIE and C-FS with high frequency. It is seen that the pupil–teacher ratio, the proportion of blacks, and the per capita crime rate are the key factors affecting housing prices. However, FAR misses the above important variables. $\left(2\right)$ For the value of SNV, both our method, GRIE, and FAR are 0, which means that they can successfully exclude all artificial variables.

In summary, compared with C-FS and FAR, our method has the smallest A-PE, the smallest SNV, and a simple model, which implies our method has better performance in feature screening under high-dimensional settings.

5.2. Arabidopsis thaliana Gene Data

We now turn to Arabidopsis thaliana gene data to illustrate the screening performance of our method. This dataset was developed by Wille et al. (2004) [23], who detected modules of closely connected isoprenoid genes in Arabidopsis thaliana. It is available on the website https://www.ncbi.nlm.nih.gov/pmc/articles/PMC545783 (accessed on 16 November 2022), which is composed of 834 genes from 58 different pathways in 118 samples. Chen et al. [24] found that GGPPS11 played an essential role in the generation of GGPP, which is the common precursor of several biologically important compounds (such as carotenoids, chlorophylls, and gibberellins), in Arabidopsis. Our goal is to identify the remaining 833 genes’ effects on the expression value of gene GGPPS11.

Followed by Wille et al. [23], the downloaded data $R=\{y,x_1,\dots ,x_{833}\}$ were converted to permille data by taking $1000R$. To get the original dataset, we model $0.001R$ here and consider the corresponding nonparametric additive models:

$$y=\sum _{j=1}^{833}{m}_j\left(x_j\right)+ϵ,$$

(9)

where y is the expression value of gene GGPPS11, and $\{x_1,\dots ,x_{833}\}$ are the expression values of the remaining 833 genes. Next, we adopt the above additive model on the full dataset to identify the important variables by the three mentioned methods. The results are given as follows:

(i)

Under FAR, we get one gene $\left\{x_{72}\right\}$, denoted by“model ($A_2$)”;

(ii)

Under GRIE, three genes $\{x_{140},x_{571},x_{560}\}$ are chosen, denoted by “model ($B_2$)”;

(iii)

Under C-FS, there nine genes were chosen, which are $\{x_{72},x_{105},x_{191},x_{476},x_{510},x_{517},$ $x_{554},x_{658},x_{800}\}$, and it is denoted by “model ($C_2$)”.

Again, using the nondegenerative Vuong test from Liao and Shi [22], we compare models ($A_2$) and ($B_2$). The corresponding p-value of the test is $0.012$, indicating that the above two models are not equivalent at the $5\%$ significance level. Then, we also compare models ($B_2$) with ($C_2$), and the p-value is 0. Hence, models ($B_2$) and ($C_2$) are not equivalent at the $5\%$ significance level.

Lastly, similarly to the first real data case, we compare FAR, C-FS, and GRIE through their prediction errors. Again, we randomly divide the full dataset into the training and validation sets with a ratio of 4:1 and repeat this process 100 times. Here, we also centralize the response variable y and set ${\kappa }_n=3$. For this real data, we consider the average numbers of model size and A-PE to evaluate the performance of the three models. The results are shown in

Table 9. Average numbers of model size, A-PE over 100 repetitions, and their robust standard deviations (in parentheses) of Arabidopsis thaliana gene data.

Approach	Model Size	A-PE
FAR	1.00 (0.00)	0.289 (0.099)
C-FS	10.15 (3.34)	0.282 (0.181)
GRIE	1.76 (1.18)	0.276 (0.093)

. Thus, we conclude that our proposed method has the smallest model size with the strongest ability for prediction and outstanding performance in identifying important covariates compared with the other two methods.

6. Conclusions

In this paper, we propose a novel variable screening screener (GRIE) for high-dimensional nonparametric additive models, which is a combination of the nonparametric smoothing ridge estimation and the group information. We note that our paper is one of the first to focus on the free marginal correlation assumption. Without the marginal correlation assumption, the proposed screener can totally separate the unimportant and important variables with a probability tending to one. Compared with iterative sure independence screening and forward screening, the proposed screener could essentially eliminate the computational burden and achieve strong, sure screening consistency. Furthermore, it allows the covariates to be strongly correlated and performs better than its alternative competitors. For these reasons, combining the strong, sure screening property of GRIE with the model selection property of EBIC, we propose the GRIE-EBIC method to further eliminate the noise variables and improve the accuracy of model selection. Theoretically, we establish the strong consistency of model selection for the GRIE-EBIC method, which reveals that our proposed method achieves the ideal model selection results.

We conclude this paper with a discussion of directions for future research. One direction to consider is nonparametric additive models with interaction effects between covariates, which are defined as

$$E(y\mid \mathbf{x})=\sum _{1\le j<k\le p_n}{m}_{j,k}(x_j,x_k),$$

where $x_j$ is the jth element of $\mathbf{x}$. They are the generalization of linear models with two-way interaction effects [25] that are more flexible for capturing the intersection between covariates. One potential approach may be to use the tensor splines bases to approximate each nonparametric function $m_{j,k}(·,·)$. The other direction is to study how to apply our methodology in the nonparametric generalized additive models [26,27]. The nonparametric generalized additive models admits

$$G\left\{E(y\mid \mathbf{x})\right\}=\sum _{j=1}^{p_n}{m}_j\left(x_j\right),$$

where $x_j$ is the jth elements of $\mathbf{x}$, and $G(·)$ is the link function. Since the nonparametric smoothing ridge estimation has outstanding performance in nonparametric additive models, its performance in generalized additive models may be worth investigating.