Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
100,926 characters · 20 sections · 59 citation commands
A One-Covariate-at-a-Time Multiple Testing Approach to Variable Selection in Additive Models
Variable selection has been playing a pivotal role in econometrics and statistics for statistical learning and scientific discoveries. Notable early contributions include Akaike1973, Akaike1974, and Schwartz1978. These authors suggested a unified approach to model selection, viz., choosing a parameter vector by minimizing the conventional criterion function plus an $L_{0}$ penalty to penalize the model size. However, these approaches are not feasible in high-dimensional settings. The seminal work of Tibshirani1996\ addressed this important issue by substituting the $L_{0}$ penalty with an $L_{1}$ penalty, and it has sparked extensive studies in both statistics and econometrics. Important contributions in the statistics literature include FanLi2001, ZhouHastie2005, Zou2006, FanLv2008, ZouLi2008, Zhang2010, FanFengSong, FanLv2013, and FanTang2013. Important contributions in the econometric literature include BelloniEtal2012, BelloniChernozhukov2013, BelloniChernozhukovHansen2014, BelloniEtal2017,\ ChernozhukovEtal2018, and Chudik_K_P. For a comprehensive review, see FanLiZhangZou2020.
In this paper, we propose a multiple testing approach to variable selection for high-dimensional nonparametric additive models. Recently Chudik_K_P (CKP hereafter) have proposed a One-Covariate-at-a-time Multiple Testing (OCMT) approach for linear regression models. CKP suggest regressing the dependent variable on each independent variable separately, retaining only those variables that exhibit a high correlation with the dependent variable. This strategy is often referred to as the \textquotedblleftscreening approach\textquotedblright. CKP's main contributions are twofold. First, they propose a criterion for variable selection by controlling the probability of choosing all the signals (or with some pseudo-signals) in the model. Second, unlike the usual screening approach that may miss some important \textquotedblleft hidden\textquotedblright\ signals whose net effects on the dependent variable are small, the CKP's OCMT procedure is able to pick up hidden signals with very high probability. The OCMT procedure has been applied in various applications; see, e.g., Kozbur2020, Chudik_P_S, and AhmendPesaran2022. In particular, Kozbur2020 considers a testing-based forward model selection (TBFMS) procedure in linear regression models that inductively selects covariates to add predictive power into a working statistical model. But this latter paper mainly focuses on the error bound and shows that the proposed procedure is able to achieve estimation rates matching those of Least Absolute Shrinkage and Selection Operator (Lasso) and post-Lasso. Furthermore, Sharifvaghefi2023 extends OCMT to cases with many highly correlated covariates and allows the number of pseudo signals to grow at the same rate as the sample size.
Our paper contributes to the literature by extending the CKP's\ OCMT approach from parametric models to nonparametric additive models. Like CKP, we estimate the net effect of each variable on the dependent variable one by one, possibly with some preselected variables. The selected variables are those whose net effects exceed some threshold value, specified to ensure the probability of selecting all the signals is very high. The statistics constructed in this paper are much more complicated than the $t$-statistics in CKP and might not even exhibit a well-behaved limiting distribution. In addition to investigating a different model, our paper differs from that of CKP in some other important aspects. First, we generalize the definition of hidden signals (in Table 2 in Section (ref)). Technical details are updated accordingly. Second, CKP chose tuning parameters similar to a Bonferroni correction of the cumulative distribution function of the standard normal. This choice implicitly requires a certain degree of approximation of the standard normal to the $t$-statistic distribution in their study. Instead, we select the tuning parameters using the classic Bayesian information criterion (BIC). Third, we add an adaptive group-Lasso-based post-OCMT step to eliminate pseudo-signals that cannot be eliminated with very high probability in CKP. This step adds very little computation burden because the dimension is reduced dramatically before the last-step estimation. This additional step is not needed in theory, but it aims at eliminating the pseudo-signals and thereby enhances the out-of-sample forecasting performance in practice.
One competing method to ours in the literature is the adaptive group Lasso (AGLASSO) proposed by HuangEtal2010. The AGLASSO adds some adaptive penalty term to the usual least squares loss function in the spirit of Zou2006. Even though we also use AGLASSO to eliminate the pseudo-signals after the OCMT procedure, our approach allows much faster computation and provides more reliable estimates than their approach.
We consider various setups in the simulation studies and compare the above post-OCMT AGLASSO procedure with that based on the OCMT\ alone or the AGLASSO procedure of HuangEtal2010 alone. We find that the former one generally outperforms the latter two significantly. We apply our method on a dataset from the Longitudinal Survey on Rural Urban Migration in China (RUMiC) and the empirical results also demonstrate the excellent performance of our procedure in finite samples.
The remainder of the paper is structured as follows. In the next section, we illustrate our approach through a single-stage procedure that is silent to hidden signals. In Section (ref), we present the more powerful multiple-stage procedure. We investigate the finite sample properties of our procedure through Monte Carlo experiments in Section (ref) and an empirical application in Section (ref). We conclude the paper in Section (ref). The proofs of all propositions and theorems in the paper are relegated to Appendix (ref). The online supplement contains some additional technical materials that include the proofs of the technical lemmas in Appendix (ref) and some additional results in the simulation and application. To facilitate reading, We present our procedure in detail for practitioners in Appendix (ref) and the idea of the proofs in Appendices (ref) and (ref).
Notation. For a generic real matrix $\boldsymbol{A=}\left\{ a_{ij}\right\} \boldsymbol{,}$ let $\left\Vert \boldsymbol{A}\right\Vert =\left[ \lambda_{\max}\left( \boldsymbol{A}^{\prime}\boldsymbol{A}\right) \right] ^{1/2}$ denote the spectral norm and $\left\Vert \boldsymbol{A} \right\Vert _{\infty}=\max_{ij}\left\vert a_{ij}\right\vert .$ When $\boldsymbol{A}$ is symmetric, $\lambda_{\max}\left( \boldsymbol{A}\right) \ $and $\lambda_{\min}\left( \boldsymbol{A}\right) \ $denote its maximum and minimum eigenvalues, respectively. For vector $\boldsymbol{x},$\ $\left\Vert \boldsymbol{x}\right\Vert $ denotes its Euclidean norm. For the deterministic series $\left\{ a_{n},b_{n}\right\} _{n=1}^{\infty}$, we denote $a_{n}\propto b_{n}$ if $0<C_{1}\leq\lim\inf_{n\rightarrow\infty}\left\vert a_{n}/b_{n}\right\vert \leq\lim\sup_{n\rightarrow\infty}\left\vert a_{n} /b_{n}\right\vert \leq C_{2}<\infty$ for some constants $C_{1}$ and $C_{2},$ $a_{n}\lesssim b_{n}$ if $\lim\sup_{n\rightarrow\infty}\left\vert a_{n} /b_{n}\right\vert \leq C<\infty$ for some $C$ that does not depend on $n,$ $a_{n}\gtrsim b_{n}$ if $b_{n}\lesssim a_{n},$ $a_{n}\ll b_{n}$ if $a_{n}=o\left( b_{n}\right) ,$ and $a_{n}\gg b_{n}$ if $b_{n}\ll a_{n}.$ $\mathcal{A}^{c}$ denotes the complement of the set $\mathcal{A}$. $\overset{P}{\rightarrow}$ denotes convergence in probability. $C$ and $M$ denote some positive constants that may vary from line to line.
In this section we give the model and definitions of various versions of signals and noises. Then we will provide the one-stage procedure for variable selection, present the basic assumptions and study the asymptotic properties of our one-stage procedure.
Recently, CKP proposed a powerful multiple-stage procedure for linear models. Our paper aims to generalize the results of CKP to the nonparametric additive models. The notations and technical details in CKP are already quite tedious, and they would be even more so in this paper. To facilitate the exposition, we start with a simple case where there are no pre-determined variables, and we conduct only one-stage multiple testing. We will present the more powerful multiple-stage procedure and show its validity in Section (ref).
Suppose the model is
where $Y$ is the dependent variable, $X_{1},$ $X_{2},\ldots,$ and $X_{p^{\ast }}$ are random independent variables, $\varepsilon$ is an unobserved error term, and $f^{\ast}$ is an unknown smooth function. Even though we have only $p^{\ast}$ signal variables, namely $X_{1},$ $X_{2},\ldots,$ $X_{p^{\ast}}$, that should be included into the regression model in ((ref)), we do not know this truth before the data reveal the fact. The realistic situation is that the $p^{\ast}$ signals are contained in a set $\mathcal{S} _{n}=\left\{ X_{j},j=1,2,\ldots,p_{n}\right\} ,$ where $p_{n}$ can be much larger than the sample size $n$ and $p_{n}\propto n^{B_{p}}$ for some $B_{p}>0$. $\mathcal{S}_{n}$ is the set for all candidate variables, and we refer to it as the active set. We assume that $E\left( \varepsilon|X_{1},X_{2},\ldots,X_{p^{\ast}},X_{p^{\ast}+1},\ldots,X_{p_{n} }\right) =0.$ The target of the model selection is to pick up those signals among $\mathcal{S}_{n}.$
To avoid the curse of dimensionality in nonparametric estimation, we impose the additive structure on $f^{\ast}$, that is,
We allow the additive components $f_{l}^{\ast}\left( X_{l}\right) $ to change with $n$ but we suppress the dependence of $f_{l}^{\ast}\left( X_{l}\right) $ on $n$ for notional convenience. Obviously, the individual functions $f_{j}^{\ast},$ $j=1,\ldots,p^{\ast},$ cannot be identified without certain suitable normalizations. We impose the normalization by assuming that $E[f_{j}^{\ast}\left( X_{j}\right) ]=0$ for each $j$ and rewrite the model as
where $\mu=E\left( Y\right) .$ Since $\mu$ can be estimated by the sample mean of the $Y$ variable at the usual $\sqrt{n}$-rate and this estimation does not affect the estimation of the nonparametric additive component, for the simplicity of presentation, we assume $\mu=0$ below.
Since we will focus on the $B$-spline-based nonparametric theory, it is standard to assume compact support for each regressor (see, e.g., Horowitz_Mammen2014, Chen_handbook and HuangEtal2010). Without loss of generality, we assume that the support for $X_{j}$ is $[0,1\mathcal{]}$ for all $j$. Let $\alpha_{1}$ be a non-negative integer, $\alpha_{2}\in(0,1],$ and $d=\alpha_{1}+\alpha_{2}.$\ Denote the class of all $\alpha_{1}$ times continuously differentiable real-valued functions on $[0,1]$ by $C^{\alpha_{1}}\left( [0,1]\right) .$\ Define $d$-th smooth real-valued functions on $[0,1]$ as
For notational simplicity, we will restrict our attention to the case where $f_{j}$'s are smooth enough and belong to $\Lambda^{d}\left( [0,1]\right) $ with $d>1$. The case of different smoothness parameters only complicates the notation but does not bring in any new insight.
The idea of one-stage procedure is that we estimate the impact of $X_{l}$ on $Y,$ $l=1,2,...,p_{n},$ one by one. So we run $p_{n}$ estimations in total and will keep those variables that are significant enough. Since in each regression we only have one covariate, we are not estimating $f_{l}^{\ast}$ when the explanatory variable is $X_{l}$. Instead, we are estimating the conditional expectation $f_{l}\left( X_{l}\right) \equiv E\left( Y|X_{l}\right) .$ We define the net impact (or the net effect) of $X_{l}$ on $Y$ as \[ \theta_{l}\equiv\left\{ E\left[ f_{l}\left( X_{l}\right) ^{2}\right] \right\} ^{1/2}=\left\{ E\left[ \left( \sum_{j=1}^{p^{\ast}}\sigma _{lj}\right) ^{2}\right] \right\} ^{1/2}, \] where $\sigma_{lj}=E[f_{j}^{\ast}\left( X_{j}\right) |X_{l}].$ Since we allow $f_{l}^{\ast}\left( X_{l}\right) $ to change with $n,$ $\theta_{l}$ might change with $n$ as well. But we suppress its dependence on $n$ for notational convenience. Here, $\sigma_{lj}$ plays the role of the scaled covariance between $X_{j}$ and $X_{l}$ for the linear model considered by CKP. To better understand the connection between the net impact defined here and in CKP, we refer the readers to the results in equations ((ref)) and ((ref)), where we approximate $f_{l}\left( X_{l}\right) $ with certain linear functions $f_{nl}\left( X_{l}\right) .$ Ignoring the bias from the approximation, one can see the connection more clearly.
Obviously, $\theta_{l}$ can be 0 or close to 0 for signals, and $\theta_{l}$ can be nonzero or large for non-signals. As in CKP, we also have four possibilities as tabulated in Table (ref). Cases (I) and (IV) in Table 1 are desirable cases. Case (III) happens when some non-signals are not independent of the signals.\ The hidden signals defined in Case (II) are rare in the linear case, and it is also rare in the nonparametric case. We generalize the definition of hidden signals to Table 2 in the next section where $\theta_{l}$ is non-zero but small relative to the sample size.
As we shall see, our one-stage procedure is silent on picking up hidden signals and eliminating pseudo-signals. For the hidden signals, we will propose a multiple-stage procedure in Section (ref)\ that can effectively pick them up. To eliminate the pseudo-signals, as a post-procedure we propose to re-estimate the model using adaptive group Lasso in Section (ref).
We assume that there are $p^{\ast\ast}$ pseudo-signals. Without loss of generality, we denote them to be \[ \left\{ X_{p^{\ast}+1},X_{p^{\ast}+2},\ldots,X_{p^{\ast}+p^{\ast\ast} }\right\} . \] Below we focus on the one-stage procedure in this section and postpone the multi-stage case to the next section.
Suppose that we have $n$ observations $\left\{ (y_{i},x_{1i},...,x_{p_{n} i})\right\} _{i=1}^{n}$ that are drawn from the distribution of $(Y,X_{1},...,X_{p_{n}}).$ The data are given in a $n\times\left( p_{n}+1\right) $ matrix \[ \left( \boldsymbol{y},\boldsymbol{x}_{1},\boldsymbol{x}_{2},\ldots ,\boldsymbol{x}_{p_{n}}\right) \] where $\boldsymbol{y}=\left( y_{1},y_{2},\ldots,y_{n}\right) ^{\prime}$ and $\boldsymbol{x}_{l}=\left( x_{l1},x_{l2},\ldots,x_{ln}\right) ^{\prime}$ for $l=1,...,p_{n}$. We propose to choose finite order (e.g., cubic) B-spline basis functions $\left\{ \psi_{l}\left( x\right) \right\} _{l=1}^{m_{n}} $\ on $\left[ 0,1\right] $\ to approximate the unknown functions $f_{l}$'s.
B-splines are piecewise-defined polynomial functions that can be used to construct curves and surfaces in numerical analysis. They offer a flexible way to model and control the shape of these curves and surfaces. A B-spline of order $n$ is a piecewise-defined polynomial function of degree $n-1$. B-splines of order one are piecewise constant functions, B-splines of order two are piecewise linear functions, B-splines of order three are piecewise quadratic functions, and so on. The points at which different polynomial pieces connect are called knots, and the knot vector specifies where these knots are. For the detailed definition and properties of B-spline bases, see Stone and deBoor. We list some properties of the B-spline basis functions in Lemma (ref). Other popular basis functions include polynomials, trigonometric polynomials, splines, and orthogonal wavelets. We refer the readers to Chen_handbook for a nice review about sieve estimation.
Since we normalize $E[f_{j}^{\ast}\left( X_{j}\right) ]=0,$ we similarly normalize the basis as \[ \phi_{jl}\left( x\right) =\psi_{j}\left( x\right) -n^{-1}\sum_{i=1} ^{n}\psi_{j}\left( x_{li}\right) . \] This is a standard practice; see, e.g., HuangEtal2010. For notational simplicity, we will write $\phi_{j}\left( x\right) $ for $\phi_{jl}\left( x\right) $. Let $P^{m_{n}}\left( x\right) =\left[ \phi_{1}\left( x\right) ,\phi_{2}\left( x\right) ,\ldots,\phi_{m_{n}}\left( x\right) \right] ^{\prime},$ an $m_{n}\times1$\ vector. Define
which are the population coefficient and the error term in the regression of $Y$ on $P^{m_{n}}\left( X_{l}\right) $. Note that we suppress the dependence of $\boldsymbol{\beta}_{l}$ on the sample size $n$.
For the one-stage procedure, we conduct the regression of $Y$ on $P^{m_{n} }\left( X_{l}\right) ,$ $l=1,2,\ldots,p_{n},$ one by one. Let $\mathbb{X} _{li}=P^{m_{n}}\left( x_{li}\right) $ be the approximating function basis at the $i$th observation for $X_{l}$.$\ $Let $\mathbb{X}_{l}=(\mathbb{X} _{l1},\mathbb{X}_{l2},\ldots,\mathbb{X}_{ln})^{\prime}$ be the $n\times m_{n}$ \textquotedblleft design\textquotedblright\ matrix for $X_{l}.$ For $X_{l},$ we regress $\boldsymbol{y}$ on $\mathbb{X}_{l}\ $to obtain \[ \boldsymbol{\hat{\beta}}_{l}=\left( \mathbb{X}_{l}^{\prime}\mathbb{X} _{l}\right) ^{-1}\mathbb{X}_{l}^{\prime}\boldsymbol{y}. \] We construct the test statistic as\footnote{As a referee has noted, one can define an alternative test statistic $\tilde{\mathcal{X}}_{l}=n\hat {\boldsymbol{\beta}}^{\prime}\hat{\boldsymbol{\beta}}$, which also works under certain rank conditions (see, e.g., Assumption (ref) below). We opted for $\hat{\mathcal{X}}_{l}$ for two reasons. First, $\hat{\mathcal{X} }_{l}$ resembles the usual chi-squared statistic under conditional homoskedasticity. Because we do not want to model the conditional heteroskedasticity of unknown form, we cannot take into account conditional heteroskedasticity explicitly in constructing the test statistic $\hat{\mathcal{X}}_{l}.$ Despite this, our asymptotic theory allows for conditional heteroskedasticity in the error term. Second and more importantly, $\hat{\mathcal{X}}_{l}$ is scale-free whereas $\tilde{\mathcal{X}}_{l}$ is not. The latter makes it very challenging to choose the range to search the constant $C$ in $\varsigma_{n}$ defined below.}
where $\hat{\sigma}_{l}^{2}=n^{-1}\sum_{i=1}^{n}\hat{u}_{li}^{2}$ and $\hat {u}_{li}$ is the residual from the above regression. Then we define the first-stage OCMT selection indictor as
where $\mathbf{1}\left( \cdot\right) $ is the usual indicator function and $\varsigma_{n}$ is a threshold value.
For the linear model in CKP with one covariate at a time, $\mathcal{\hat{X} }_{l}$ is asymptotically $\chi^{2}\left( 1\right) $ under conditional homoskedasticity, and one can follow their lead to consider threshold values for the associated $t$-statistics based on the adjusted normal critical values. Nevertheless, such a result is not available in our framework due to the divergent dimension of regressors in the sieve estimation. In addition, the potential presence of conditional heteroskedasticity greatly complicates our analysis too. What we really need is to show that $\mathcal{\hat{X}}_{l}$ behaves distinctly for signals and noises so that a suitable choice of the threshold value $\varsigma_{n}$ can help us separate the signals from the noises. For these reasons, we do not associate $\mathcal{\hat{X}}_{l}$ with any asymptotic distribution. Instead, we will set $\varsigma_{n}\propto \kappa_{n}\log\left( m_{n}\right) m_{n}$ for a positive series $\kappa_{n}$ that diverges to infinity slowly as in Assumption (ref). For more details, see the remark on Assumption (ref) in the next subsection.
To study the asymptotic properties of the one-stage procedure, we impose the following assumptions.
Assumption \ref*{A:iid} imposes an i.i.d. condition on the observations and a standard conditional moment restriction. The extension to weakly dependent observations is possible but left for future research. The generalization to independently non-identically distributed (i.n.i.d.) case is straightforward because the main inequalities in Lemmas A.1 and A.2 allow for i.n.i.d. observations. We keep using the i.i.d. assumption for notational convenience. In Assumption \ref*{A:p}, we assume that the number of signals is fixed; the number of pseudo-signals is allowed to increase as $n$ increases but at a slower rate than that of the total candidate variables. It is also possible to allow $p^{\ast}$ to diverge to infinity (see Section (ref)). Assumption \ref*{A:supp} restricts the support of $X$ to be $\left[ 0,1\right] .$ This is a very common condition for nonparametric additive models; see, e.g., Li2000 and Horowitz_Mammen2014. Assumption \ref*{A:epsilon} imposes some tail conditions $\varepsilon,$ which is also assumed in CKP but weaker than the commonly used sub-exponential condition in the variable selection literature and the one for the adaptive group Lasso in HuangEtal2010. The primary use of this condition is to derive probability bound for errors.
Assumption \ref*{A:fl} imposes fairly weak smooth condition on $f_{l}\left( x\right) $, which is weaker than the commonly used condition $d\geq2$. Note that $d>1$ is needed for Assumption \ref*{A:mn}. Assumption \ref*{A:tech} is a technical assumption needed to simplify the proof.\ Specifically, the boundedness of $f_{j}^{\ast}$ implies $U_{l}$ defined in equation ((ref)) has the same tail behavior as $\varepsilon$. The bounded conditional second moment of $\varepsilon$\ is to ensure some nice properties of $U_{l}\phi_{j}\left( X_{l}\right) $ and this assumption is also common in the sieve literature (see, e.g., Newey1997).\ This assumption is mild given that we assume all $X_{l}$'s have compact support and the tails of $\varepsilon$ decay exponentially fast. Assumption \ref*{A:mn} imposes conditions on $m_{n}.$ First, we need $B_{m}<1/3$ such that $nm_{n} ^{-3}\rightarrow\infty.$ The last condition is necessary for $\left\Vert n^{-1}\mathbb{X}_{l}^{\prime}\mathbb{X}_{l}\right\Vert \propto m_{n}^{-1}$ to hold with very high probability. To see why, note that Lemma (ref) in Appendix (ref) suggests that $\left\Vert E\left[ P^{m_{n} }\left( X_{l}\right) P^{m_{n}}\left( X_{l}\right) ^{\prime}\right] \right\Vert \propto m_{n}^{-1}.$ We need $m_{n}^{-1}\gg\left( m_{n}/n\right) ^{1/2},$ or equivalently, $nm_{n}^{-3}\rightarrow\infty,$ in order to ensure that $n^{-1}\mathbb{X}_{l}^{\prime}\mathbb{X}_{l}$ is close to $E\left[ P^{m_{n}}\left( X_{l}\right) P^{m_{n}}\left( X_{l}\right) ^{\prime }\right] $ with very high probability. Second, we need $B_{m}>1/\left( 1+2d\right) $ to ensure that the approximation bias is asymptotically negligible in comparison with the asymptotic variance term: $m_{n}^{-d} \ll\left( m_{n}/n\right) ^{1/2}$.
Assumption \ref*{A:xi_n} imposes conditions on the threshold value $\varsigma_{n}$ that ensures the separability of the signals from noises. If the true value of $\boldsymbol{\beta}_{l}$ is $\boldsymbol{0}$, $\mathcal{\hat {X}}_{l}$ in equation ((ref)) is $O_{P}\left( m_{n}\right) .$ In this case, to ensure $\widehat{\mathcal{J}}_{l}=0$ with very high probability, we can take $\varsigma_{n}\propto\kappa_{n}\log\left( m_{n}\right) m_{n}$ and lose some power up to $\kappa_{n}\log\left( m_{n}\right) .$ Here, $\kappa_{n}$ can be any series diverging to infinity slowly, e.g., $\left[ \log\left( m_{n}\right) \right] ^{\epsilon}$ for some small $\epsilon>0.$ The loss of the power to some degree is inevitable because of the nature of the multiple testing procedure when the number of tests goes to infinity. In contrast, CKP choose their threshold by the Bonferroni correction of the standard normal. We choose not to follow them because of the following reasons. First, given the divergence of $m_{n},$ $\mathcal{\hat{X}}_{l}$ does not converge to a chi-square distribution asymptotically even in the homoskedastic case so that we cannot use chi-square distribution to approximate the finite sample distribution of $\mathcal{\hat {X}}_{l}.$ Under conditional heteroskedasticity, $\mathcal{\hat{X}}_{l}$ does not converge to a chi-square distribution even if $m_{n}$ is held fixed. So our procedure does not rely on the chi-square approximation. Second, even if we can do the approximation, the cumulative density function (CDF) of a chi-square distribution is very complicated. We do not have a rate for the inverse of its CDF evaluated at a certain rate (e.g., $n^{-C}$) like the case of normal CDF. For these reasons, we do our selection based on the asymptotic results. Specifically, we will take $\varsigma_{n}=C\kappa_{n}\log\left( m_{n}\right) m_{n}$ for some $\kappa_{n}$ and choose the value of $C$ by the classic BIC. The details can be found in Appendix (ref). This $\varsigma_{n}$ helps separate the signals from the noises with very high probability, as demonstrated in the next section.\ Recall that in Assumption \ref*{A:p} we assume that $p_{n}$ go to infinity at a polynomial rate of $n,$ same as $m_{n}.$ This ensures that $\log\left( m_{n}\right) \propto\log p_{n}.$ Therefore, theoretically we only need to put $\log\left( m_{n}\right) $ in $\varsigma_{n}$ to have a control of $p^{\ast\ast}$ and $p_{n}$ for the TPR, FDR and FPR defined after Proposition (ref) below. We postpone the discussion on the comparison of technical conditions required for our procedure and those required for the AGLASSO to Section 3.2.
We present the first theoretical result in this paper. It derives the probability bounds for the \textquotedblleft Type-I\textquotedblright\ and \textquotedblleft Type-II\textquotedblright\ errors.
The proof of Proposition (ref) is tedious. We provide some technical details in Appendix (ref) before we formally prove the proposition in Appendix (ref). An implication of Proposition (ref) is that for the well-chosen threshold value $\varsigma_{n},$ the above one-stage procedure can separate the signals with $\theta_{l}\gtrsim\kappa_{n}\log\left( m_{n}\right) ^{1/2}\left( m_{n}/n\right) ^{1/2}$ from the noises with $\theta_{l}=0.$ Of course, we may have some intermediate case where $0<\theta_{l}\lesssim\log\left( m_{n}\right) ^{1/2}\left( m_{n}/n\right) ^{1/2}$ for which the above procedure fails to do so. Note that the convergence rate for each additive term in Stone is $\left( m_{n}/n\right) ^{1/2}.$ Therefore, our procedure loses some power up to the order of $\log\left( n\right) $, a common scenario in the variable selection literature.
Following the literature and CKP, we define the true positive rates (TPR) and the false positive rates (FPR) respectively as
Based on the test statistic defined in equation ((ref)) and its property developed in Proposition (ref), we introduce the generalized definitions of signals and noises in Table (ref). With this definition of hidden signals, our result in Section (ref) provides theoretical justification for the necessity of a multi-stage procedure capable of detecting signals not identified in the first stage, yet with non-zero net effects.
With a little abuse of notation, we continue to use $p^{\ast}$ and $p^{\ast\ast}$ to denote the number of signals and pseudo-signals at the sample level, and assume that they satisfy Assumption \ref*{A:p}. Adopting the definitions in Table 2, we define the false discovery rates (FDR) as \[ \text{FDR}_{n}=\frac{\sum_{l=1}^{p_{n}}\mathbf{1}\left( \widehat{\mathcal{J} }_{l}=1\text{, }\left\{ E\left[ f_{l}^{\ast}\left( X_{l}\right) ^{2}\right] \right\} ^{1/2}=0,\text{ and }\theta_{l}\lesssim\log\left( m_{n}\right) ^{1/2}\left( m_{n}/n\right) ^{1/2}\right) }{\sum_{l=1} ^{p_{n}}\widehat{\mathcal{J}}_{l}+1}. \]
Apparently, the definitions in Table 2 generalize the definitions in Table 1.
The one-stage procedure is valid in terms of TPR, only if the net effects of all signals are strong enough. Specifically, we need the following assumption.
We note the definition of no hidden signals in the above assumption is equivalent to the one in Table 2, given the way $\kappa_{n}$ is defined in Assumption 8. If Assumption \ref*{A;no_hidden} fails to hold, we need the multiple-stage procedure to pick up the hidden signals.
Based on the results in Proposition (ref), we present the results for TPR$_{n}$, FPR$_{n}$, and FDR$_{n}$ in the following theorem.
Theorem (ref) implies that all of TPR$_{n},$ FPR$_{n}$ and FDR$_{n}$ can be well controlled provided we assume away hidden signals at the sample level. Note that Theorem (ref)(i)--(ii) focuses on the asymptotic properties of TPR$_{n}$ and FPR$_{n}$ while the last part of Theorem (ref) reveals that the false discovery rate is asymptotically vanishing in large samples.
In the next section, we turn to the multiple-stage procedure that does not rely on Assumption \ref*{A;no_hidden}.
In this section we propose a multiple-stage procedure to select variables for the nonparametric additive models.
As mentioned above, we may not identify a signal $X_{l}$ whose net effect satisfies $\theta_{l}\lesssim\kappa_{n}\log\left( m_{n}\right) ^{1/2}\left( m_{n}/n\right) ^{1/2},$ even in the case where the marginal effect of $X_{l}$ on $Y,$ namely, $\left\{ E\left[ f_{l}^{\ast}\left( X_{l}\right) ^{2}\right] \right\} ^{1/2},$ is large enough.\footnote{The marginal effect defined here is slightly different from that in the econometrics literature. For example, the marginal effect of $X_{l}$ on $Y$ is defined as $\beta_{l}$ in the CKP's linear model: $Y=\beta_{0}+\sum _{l=1}^{p^{\ast}}X_{l}\beta_{l}+\varepsilon,$ but it refers to $X_{l}\beta _{l}$.in this paper.} Consequently, Theorem (ref) does not hold without Assumption \ref*{A;no_hidden} which assumes away small sample hidden signals. In contrast, Lemma (ref) in Appendix (ref) underpins the result that as long as $\left\{ E\left[ f_{l}^{\ast}\left( X_{l}\right) ^{2}\right] \right\} ^{1/2}\gtrsim\kappa_{n}\log\left( m_{n}\right) ^{1/2}\left( m_{n}/n\right) ^{1/2}$ for some slowly divergent series $\kappa_{n}$ and some full rank condition holds, the net effect of $X_{l}$ will be strong enough to be picked up at certain stage of a multiple-stage procedure. A by-product of this Lemma is the existence of at least one signal for Stage 1 (see the second part of this Lemma).
To introduce the multiple-stage procedure, we need some extra notations. Suppose at certain stage after stage 1, we have pre-selected $\iota_{n}$ variables from the active set $\mathcal{S}_{n}=\left\{ X_{j},\text{ }j=1,...,p_{n}\right\} $ based on some selection procedure to be described below. To avoid confusion, we denote these pre-selected variables as $Z_{1},\ldots,Z_{\iota_{n}}$. Let $\boldsymbol{Z}\equiv\left( Z_{1} ,\ldots,Z_{\iota_{n}}\right) ^{\prime}\ $and $P^{m_{n}}\left( \boldsymbol{Z} \right) =(P^{m_{n}}\left( Z_{1}\right) ^{\prime},\ldots,$ $P^{m_{n}}\left( Z_{\iota_{n}}\right) ^{\prime})^{\prime}.$ Note that $P^{m_{n}}\left( \boldsymbol{Z}\right) $ is an $\iota_{n}m_{n}\times1$ vector for $\boldsymbol{Z}.$ At the next stage, we consider the nonparametric additive regression of $Y$ on $\boldsymbol{Z}$ and an $X_{l}$ that has not been selected so far and we do this one by one for all $p_{n}-\iota_{n}$ non-selected variables $X_{l}.$ We define the impact of $X_{l}$ on $Y$ after controlling $\boldsymbol{Z}$ as
where $\Phi_{\boldsymbol{Z}}\equiv E[P^{m_{n}}\left( \boldsymbol{Z}\right) P^{m_{n}}\left( \boldsymbol{Z}\right) ^{\prime}]$ and $\mu _{lj,\boldsymbol{Z}}\equiv E\left\{ \left. f_{j}^{\ast}\left( X_{j}\right) -P^{m_{n}}\left( \boldsymbol{Z}\right) ^{\prime}\Phi_{\boldsymbol{Z}} ^{-1}E\left[ P^{m_{n}}\left( \boldsymbol{Z}\right) f_{j}^{\ast}\left( X_{j}\right) \right] \right\vert X_{l}\right\} $ denotes the effect of $X_{l}$ on $f_{j}^{\ast}\left( X_{j}\right) $ after controlling the effects of $\boldsymbol{Z}$. Apparently, we suppress the dependence of $\theta _{l,\boldsymbol{Z}}$ on the sample size $n$.
At the sample level, let $\mathbb{Z}$ and $\mathbb{X}_{l}$ denote the $n\times m_{n}\iota_{n}$ and $n\times m_{n}$ \textquotedblleft design matrices\textquotedblright\ for $\boldsymbol{Z}$ and $X_{l},$ respectively. That is,
where $\mathbb{Z}_{l}=\left( \mathbb{Z}_{l1},\ldots,\mathbb{Z}_{ln}\right) ^{\prime}\ $is a $n\times m_{n}$ matrix, $\mathbb{Z}_{li}=P^{m_{n}}\left( z_{li}\right) \ $and $\mathbb{X}_{li}=P^{m_{n}}\left( x_{li}\right) .$ Define $M_{\mathbb{Z}}\equiv I_{n}-\mathbb{Z}\left( \mathbb{Z}^{\prime }\mathbb{Z}\right) ^{-1}\mathbb{Z}^{\prime}.$ By the result of partitioned regressions, the coefficient of $P^{m_{n}}\left( X_{l}\right) \ $is estimated by \[ \boldsymbol{\hat{\beta}}_{l,\boldsymbol{Z}}=\left( \mathbb{X}_{l}^{\prime }M_{\mathbb{Z}}\mathbb{X}_{l}\right) ^{-1}\mathbb{X}_{l}^{\prime }M_{\mathbb{Z}}\boldsymbol{y}. \] To determine whether $X_{l}$ should be treated as a signal variable, we propose the following test statistic
where $\hat{\sigma}_{l,\boldsymbol{Z}}^{2}=n^{-1}\sum_{i=1}^{n}\hat {\varepsilon}_{li}^{2}$ and $\hat{\varepsilon}_{li}$ is the residual from the regression $\boldsymbol{y}$ on $\left( \mathbb{Z},\mathbb{X}_{l}\right) $.
We will study the asymptotic properties of $\mathcal{\hat{X}} _{l,\boldsymbol{Z}}$ in the next subsection which lay down the foundation for our multiple-stage procedure.
To study the asymptotic properties of $\mathcal{\hat{X}}_{l,\boldsymbol{Z}},$ we impose the following technical conditions.
\noindentAssumption \ref*{A:p}' $p^{\ast}\ $is a positive integer that does not vary with $n.$\ $p^{\ast\ast}\lesssim n^{B_{p^{\ast\ast }}}\ $and $p_{n}\propto n^{B_{p}}$ for some $B_{p}>B_{p^{\ast\ast}}\geq0.$ Further, $B_{p^{\ast\ast}}<\left( 1-3B_{m}\right) /2.$
\noindentAssumption \ref*{A:fl}' $f_{l}\left( \cdot\right) =E\left( Y|X_{l}=\cdot\right) \in\Lambda^{d}\left( [0,1]\right) \ $with $d>1$ for $l=1,\ldots,p_{n}$. $f_{j}^{\ast}\in\Lambda^{d}\left( [0,1]\right) \ $with $d>1$ for $j=1,\ldots,p^{\ast}$.
\noindentAssumption \ref*{A;no_hidden}' $\left\{ E[f_{j}^{\ast}\left( X_{j}\right) ^{2}]\right\} ^{1/2}\gtrsim\kappa_{n} \log\left( m_{n}\right) ^{1/2}\left( m_{n}/n\right) ^{1/2}$ for some slowly divergent series $\kappa_{n}$ as in Assumption (ref)\ and for $j=1,\ldots,p^{\ast}$.
Assumption \ref*{A:p}' strengthens Assumption \ref*{A:p} by adding one more condition on $B_{p^{\ast\ast}}$. It is imposed to ensure the good property of $\boldsymbol{\tilde{u}}_{l,\boldsymbol{Z}}$ defined in equation ((ref)) (see Lemma (ref)). This condition can be very restrictive on $B_{p^{\ast\ast}}.$ For example, when $B_{m}=1/4,$ this condition implies $B_{p^{\ast\ast}}<1/8$. Assumption \ref*{A:fl}' strengthens Assumption \ref*{A:fl} so that equation ((ref)) in Appendix A.2 holds. As remarked at the beginning of last subsection, we do not impose Assumption \ref*{A;no_hidden} for the multiple-stage procedure, as long as the marginal effect of the signal is not too weak as imposed in Assumption \ref*{A;no_hidden}'. This assumption seems inevitable. It is analogous to Assumption 6 in CKP for the linear regression model and similar to the so-called `beta-min' condition that is commonly assumed in the penalized regression literature (see, e.g., Chapter 7.4 of Buhlmann2011). Assumption \ref*{A:full_rank2} is the common rank condition for nonparametric additive regressions. We also require it to hold for the case when we add one noise variable. This seems inevitable because we will run the regression with regressors being all signals and pseudo-signals plus one noise variable in the multiple-stage procedure with very high probability.\ Assumption \ref*{A:noisevariable} is also inevitable and implicitly imposed in CKP for the linear regression models.
It is worth mentioning that Assumption \ref*{A:full_rank2} plays a similar role to the \textquotedblleft restrictive eigenvalues\textquotedblright \ condition in BickelRitovTsyvakov2009 and BelloniEtal2012. The restrictive-eigenvalues\ condition requires certain full rank conditions on all possible design matrices composed of a certain number of covariates. In contrast, our OCMT only requires full rank conditions on the design matrices composed of signals and pseudo-signals, and permits arbitrary correlations among the noise variables. Obviously, these two sets of conditions are non-nested. In addition, HuangEtal2010 impose essentially the same set of assumptions except the rank conditions discussed here. The main advantage of the OCMT is that it does not require any numerical min-search of an objective function, can be computed much faster, and deliver more reliable results.
The following proposition presents the probability bounds for the \textquotedblleft Type-I\textquotedblright\ and \textquotedblleft Type-II\textquotedblright\ errors when we have some pre-selected variables.
The proof of Proposition (ref) is rather tedious. We provide some technical discussions on the proof in Appendix (ref) before we formally prove it in Appendix (ref).
Like Proposition (ref), Proposition (ref) implies that for the well-chosen threshold value $\varsigma_{n},$ the use of the test statistic $\mathcal{\hat{X}}_{l,\boldsymbol{Z}}$ helps to separate variables with large value of $\theta_{l,\boldsymbol{Z}}$ from those with small value of $\theta_{l,\boldsymbol{Z}}.$ This observation will be used in our multiple-stage procedure to select all signal variables.
We present the multiple-stage procedure as follows.
We conduct the first-stage selection as in Section (ref) by constructing the test statistic $\mathcal{\hat{X}}_{l}$ as in equation ((ref)) and using the threshold value $\varsigma_{n}$ that satisfies the condition in Assumption \ref*{A:xi_n}. We re-label the selection indicator in equation ((ref)) as \[ \widehat{\mathcal{J}}_{l,\left( 1\right) }=\mathbf{1}\left( \mathcal{\hat {X}}_{l}>\varsigma_{n}\right) \text{ for }l=1,2,\ldots,p_{n}. \] We collect all the variables selected in stage 1 into the vector $\boldsymbol{Z}_{\left( 1\right) }$, and denote the index set of the selected variables by $S_{\left( 1\right) }.$ For the second stage, we denote the index set of the active variables in stage 2 by $\Psi_{\left( 2\right) }$ where $\Psi_{\left( 2\right) }=\left\{ 1,2,\ldots ,p_{n}\right\} \backslash S_{\left( 1\right) }$. In the second stage, we regress $Y$ on $P^{m_{n}}\left( X_{l}\right) $ with $P^{m_{n}}\left( \boldsymbol{Z}_{\left( 1\right) }\right) $ as pre-selected variables one by one for $l\in\Psi_{\left( 2\right) }.$ We construct the test statistic $\mathcal{\hat{X}}_{l,\boldsymbol{Z}_{\left( 1\right) }}$ as in equation ((ref)). We select the variable $X_{l}$ if $\widehat{\mathcal{J} }_{l,\left( 2\right) }=1,$ where \[ \widehat{\mathcal{J}}_{l,\left( 2\right) }=\mathbf{1}\left( \mathcal{\hat {X}}_{l,\boldsymbol{Z}_{\left( 1\right) }}>\varsigma_{n}\right) \text{ for }l\in\Psi_{\left( 2\right) }. \] We add all the variables selected in stage 2 into the set of variables selected in stage 1 as a new vector, and we denote it as $\boldsymbol{Z} _{\left( 2\right) }.$ We denote the index set of the selected variables ($\boldsymbol{Z}_{\left( 2\right) }$) by $S_{\left( 2\right) }$ and the index set of the active variables for stage 3 as $\Psi_{\left( 3\right) },$ where $\Psi_{\left( 3\right) }=\left\{ 1,2,\ldots,p_{n}\right\} \backslash S_{\left( 2\right) }.$ And so on and so forth. For stage $k,$ we denote the pre-selected variables as $\boldsymbol{Z}_{\left( k-1\right) },$ and the index set of the active variables as $\Psi_{\left( k\right) }.$ Then we regress $Y$ on $P^{m_{n}}\left( X_{l}\right) $ with $P^{m_{n}}\left( \boldsymbol{Z}_{\left( k-1\right) }\right) $ as pre-selected variables one by one for $l\in\Psi_{\left( k\right) }.$ We construct the test statistic $\mathcal{\hat{X}}_{l,\boldsymbol{Z}_{\left( k-1\right) }}$ as in equation ((ref)). We select the variable $X_{l}$ if $\widehat{\mathcal{J} }_{l,\left( k\right) }=1,$ where \[ \widehat{\mathcal{J}}_{l,\left( k\right) }=\mathbf{1}\left( \mathcal{\hat {X}}_{l,\boldsymbol{Z}_{\left( k-1\right) }}>\varsigma_{n}\right) \text{ for }l\in\Psi_{\left( k\right) }. \] We add all the variables selected in stage $k$ into the set of variables selected in stage $k-1$ as a new vector, and we denote it as $\boldsymbol{Z} _{\left( k\right) }.$ We stop the procedure at a stage in which no new variables are selected. We denote the stage, in which one or more variables are selected but no new variables are selected after that, as $\hat{k}_{s}$. So the OCMT procedures stops after stage $\hat{k}_{s}.$ The selection indicator for variable $X_{l}$ of the OCMT procedure is defined as follows
By construction, $\widehat{\mathcal{J}}_{l}$ is either 1 or 0. It takes value $1$ if $X_{l}$ is selected in the OCMT\ procedure and 0 otherwise.
The following theorem mainly studies the asymptotic properties of the multiple-stage procedure in terms of TPR, FPR\ and FDR.
Theorem (ref)(i) implies the OCMT procedure can terminate at step $p^{\ast}$ with very high probability. Theorem (ref) (ii)-(iv) implies that all of TPR$_{n},$ FPR$_{n}$ and FDR$_{n}$ can be well controlled. Of course, when $n\rightarrow\infty,$ we need to conduct at most $p^{\ast}$-stage procedure to determine all signals and eliminate all noise variables, with very high probability.
Because of the nature of the OCMT procedure, pseudo-signals cannot be excluded from the selection list with high probability. To eliminate pseudo-signals, we propose to employ the adaptive group Lasso after the OCMT. Then by the properties of the adaptive group Lasso, the event that the pseudo-signals are excluded and the signals are kept occurs with probability approaching one (w.p.a.1). We present this result in Theorem (ref) below. We denote the set of the variables selected by the OCMT procedure as $\hat {S}_{\text{OCMT}}\equiv S_{\left( \hat{k}_{s}\right) },$ collect them into a vector $\boldsymbol{Z}_{\text{OCMT}},$ and denote its dimension as $\hat {p}_{\text{OCMT}}$. Similarly, let $\mathbb{Z}_{\text{OCMT}}$ ($n\times\hat {p}_{\text{OCMT}}m_{n}$) denote the design matrix for $\boldsymbol{Z} _{\text{OCMT}}$ as in equation ((ref)). The post-OCMT adaptive group Lasso procedure goes as follows:
The post-selection estimation proceeds as follows. Denote the selected regressor from the above procedure as $\boldsymbol{Z}_{\text{AGLASSO}},$ and similarly denote its B-spline basis and design matrix as $P^{m_{n}}\left( \boldsymbol{Z}_{\text{AGLASSO}}\right) $ and $\mathbb{Z}_{\text{AGLASSO}},$ respectively$.$ The post selection estimator is the OLS estimator of regressing $\boldsymbol{y}$ on $\mathbb{Z}_{\text{AGLASSO}}$, which is \[ \boldsymbol{\hat{\beta}}_{\text{post}}=\left( \mathbb{Z}_{\text{AGLASSO} }^{\prime}\mathbb{Z}_{\text{AGLASSO}}\right) ^{-1}\mathbb{Z}_{\text{AGLASSO} }^{\prime}\boldsymbol{y.} \] The final fitted model is \[ P^{m_{n}}\left( \boldsymbol{Z}_{\text{AGLASSO}}\right) ^{\prime }\boldsymbol{\hat{\beta}}_{\text{post}}. \]
The above theorem is almost the same as that in HuangEtal2010 including the requirements on $\lambda_{n1}$ and $\lambda_{n2}$, with the exception that the procedure starts with the covariates post OCMT. Consequently, we only need to show that the adaptive group Lasso procedure is still valid in the post-OCMT situation, and the rest follows immediately from HuangEtal2010. In particular, under Assumption \ref*{A:mn}, the biases of the post OCMT estimators of the nonparametric additive components are asymptotically negligible so that their mean square errors (MSEs) are dominated by their asymptotic variances that are of order $O\left( m_{n}/n\right) ,$ which explains the result in Theorem (ref)(ii).
Our procedure enjoys the same theoretical property as the adaptive group Lasso. After the OCMT, the dimension of candidate variables is reduced dramatically. Thus the additional computation burden applying the post-OCMT adaptive group Lasso can be almost ignored. We note that the main advantage of our procedure is fast and reliable computation, which delivers better small sample performance as shown from our simulation studies. An implication of the above theorem is that the post-OCMT adaptive group procedure improves over the post-selection estimation results in CKP (Theorem 2 in CKP) because pseudo-signals are now eliminated with very high probability.
Allowing a diverging $p^{\ast}$ (number of true signals) is possible. The only additional technical condition apart from the restriction on the speed of $p^{\ast}$ is that $\sum_{j=1}^{p^{\ast}}f_{j}^{\ast}\left( X_{j}\right) $ is uniformly bounded. Note that this condition naturally holds for a fixed $p^{\ast}$ due to the boundedness of $f_{j}^{\ast}$. The main reason for the requirement of this condition is technical: the uniform boundedness of $\sum_{j=1}^{p^{\ast}}f_{j}^{\ast}\left( X_{j}\right) $ ensures that $U_{l} $\ defined in equation ((ref)) also satisfies the exponential decayed tail condition,\footnote{For details, see the proof of Lemma (ref).} and the tail property is necessary to apply the main inequalities to obtain the probability bounds. The uniform boundedness of $\sum_{j=1}^{p^{\ast}}f_{j}^{\ast}\left( X_{j}\right) $ was also imposed in FanFengSong$.$ Propositions (ref) and (ref) are on individual $X_{l},$\ but we do need Assumption \ref*{A:p}\textquotedblright \ on $p^{\ast}$ so that Proposition (ref) holds. It is to ensure that we can have precise estimation on a diverging design matrix.
\noindentAssumption \ref*{A:p}\textquotedblright $p^{\ast }\lesssim n^{B_{p^{\ast}}}$ for some $B_{p^{\ast}}\geq0.$\ $p^{\ast\ast }\lesssim n^{B_{p^{\ast\ast}}}\ $and $p_{n}\propto n^{B_{p}}$ for some $B_{p}>B_{p^{\ast\ast}},B_{p^{\ast}}\geq0.$ Further, $B_{p^{\ast}} +B_{p^{\ast\ast}}<\left( 1-3B_{m}\right) /2.$
We present the main results in the following theorem.
In the next section, we investigate the small sample performance of our procedure by means of Monte Carlo experiment.
To investigate the finite-sample performance of our procedure, we conduct Monte Carlo experiments in this section.
Following HuangEtal2010, we consider the following data generating processes (DGPs). In what follows, we assume
For the errors, we assume $\varepsilon\sim$ i.i.d. $N\left( 0,1\right) $ for DGPs 1--6 and 9--10. In DGPs 7--8, we check the impact of heteroskedastic errors on our methods. Specifically, we add a heteroskedastic error to the simplest and the most complicated designs in DGPs 1--6 to form DGPs 7 and 8, respectively. In DGPs 9--10, we consider the case with additive components of binary variables that mimic the application in Section (ref).
DGP 1: Four independent signals only. $Y$ is generated as follows:
Note the coefficients before $f_{i}$'s\ are set to make each signal have the same strength in terms of variance\ for independent uniform $X_{1} ,\ldots,X_{4}$. The covariates are generated as follows$:$ \[ X_{j}=W_{j}\text{ for }j=1,\ldots,4,\text{ and }X_{j}=\frac{W_{j}+U_{1}} {2}\text{for }j\geq5, \] where $W_{j},$ $j=1,\ldots,p_{n}$, and $U_{1}$ are all independent draws from $U(0,1)$. Thus, $p^{\ast}=4$ and $p^{\ast\ast}=0$ for DGP 1.\ Define the Signal-to-noise ratio to be $r_{sn}=\frac{\mathtt{sd}\left( f\right) }{\mathtt{sd}\left( \varepsilon\right) }$, and $r_{sn}=1.5$ for DGP 1.
DGP 2: Four independent signals and two pseudo-signals. $Y$ is generated from equation ((ref)). The covariates are generated as follows$:$
where $W_{j},$ $j=1,\ldots,p_{n}-2$, $U_{1},$ and $U_{2}$\ are all independent draws from $U(0,1)$. Thus, $p^{\ast}=4$ and $p^{\ast\ast}=2$ for DGP 2.
DGP 3: Four signals, and one hidden signal. $Y$ is generated from
where $f_{5}\left( X_{5}\right) =-\mathbb{E}\left[ 2.55f_{1}\left( X_{1}\right) +2.57f_{2}\left( X_{2}\right) +1.68f_{3}\left( X_{3}\right) +f_{4}\left( X_{4}\right) |X_{5}\right] .$ The covariates are generated as follows $:$
where $W_{j},$ $j=1,\ldots,p_{n}-1,$ $U_{1},$ and $U_{2}$\ are independent draws from $U(0,1)$. Then, $p^{\ast}=5$ and $p^{\ast\ast}=0$ for DGP 3, and the fifth signal is hidden by our definition.\footnote{By the distribution of covariates,
}
DGP 4: Four signals, two pseudo-signals, and one hidden signal. $Y$ is generated from equation ((ref)). The covariates are generated as follows $:$
where $W_{j},$ $j=1,\ldots,p_{n}-3,$ $U_{1},$ $U_{2},$ and $U_{3}$\ are independent draws from $U(0,1)$. Then, $p^{\ast}=5$ and $p^{\ast\ast}=2$ for DGP 4, and the fifth signal is a hidden signal.
DGP 5: Four correlated signals. $Y$ is generated from equation ((ref)). The covariates are generated as follows$:$
where $W_{j},$ $j=1,\ldots,p_{n}$, $U_{1},$ and $U_{2}$\ are independent draws from $U(0,1)$. Thus, four signals are correlated with each other, and $p^{\ast}=4$ and $p^{\ast\ast}=0$ for DGP 5.
DGP 6: Four signals, many pseudo-signals, and one hidden signal. $Y$ is generated from equation ((ref)). The covariates are generated as follows $:$
where $W_{j},$ $j=1,\ldots,p_{n}-1,$ and $U_{1}$\ are independent draws from $U(0,1)$. Then, $p^{\ast}=5$ with one hidden signal for DGP 6.
DGP 7: Four independent signals with heteroskedastic errors. $Y$ is generated from equation ((ref)) with the same covariates as in DGP 1. We assume that conditioning on $X,$ $\varepsilon$ is normal with mean 0 and variance $0.436[1+\left( X_{1}+X_{2}+X_{3} +X_{4}\right) /4]^{2},$ and the unconditional variance of $\varepsilon$ is approximately 1.
DGP 8: Four signals, many pseudo-signals, and one hidden signal with heteroskedastic errors. $Y$ is generated from equation ((ref) ) with the same covariates as in DGP 6. We assume that conditioning on $X,$ $\varepsilon$ is normal with mean 0 and variance $0.436[1+\left( X_{1} +X_{2}+X_{3}+X_{4}\right) /4]^{2}$.
DGP 9: Four independent signals with some binary variables. To mimic the application, we consider the situation with some binary covariates. Note that any function of a binary covariate can at most take two values. Without loss of generality, we focus on the case in which those binary covariates enter the model linearly. It is easy to see that our theoretical results continue to hold in the presence of some linear additive components with the main difference that they do not exhibit any approximation bias.\footnote{In this case, the test statistics for the additive linear terms are the squares of t-statistics, and the inequality continues to hold by setting, for example, $\varsigma_{n}\propto\left[ \log p_{n}\right] ^{1.1}$ for the case when $p_{n}$ is a polynomial of $n.$} We set the threshold as $\varsigma _{n}=C\left[ \log p_{n}\right] ^{1.1}$. $Y$ is generated from \[ Y=2.57f_{2}\left( X_{1}\right) +1.68f_{3}\left( X_{2}\right) +1.47X_{3}+1.47X_{4}+\varepsilon, \] where $X_{j}=W_{j},$ $j=1,2,$ $X_{j}=V_{j-2},$ $j=3,4,$ $W_{1}$ and $W_{2}$ are independent $U(0,1),$ and $V_{1}$ and $V_{2}$ are independent Bernoulli random variables with equal chances of taking value 0 or 1. The remaining covariates are generated as follows: \[ X_{j}=\frac{W_{j-2}+U_{1}}{2}\text{ for }j=5,6,...,\frac{p_{n}}{2},\text{ and }X_{j}=V_{j-\frac{p_{n}}{2}+2}\text{ for }j=\frac{p_{n}}{2}+1,...,p_{n}, \] where $W_{j},$ $j=5,6,...,\frac{p_{n}}{2}-2,$ and $U_{1},$ are independent $U(0,1),$ $V_{j},$ $j=3,...,\frac{p_{n}}{2}+2$ are distributed the same as $V_{1}.$ All $W$s$,V$s, and $U$ are independent of each other.
DGP 10: Four signals with one hidden signal in the presence of some binary variables. For this DGP, $Y$ is generated from \[ Y=2.57f_{2}\left( X_{1}\right) +1.5X_{2}+1.5X_{3}-X_{4}+\varepsilon, \] where $X_{1}=W_{1}$ is a $U(0,1),$ $X_{j}=V_{j-1},$ $j=2,3,4,$ are Bernoulli random variables with equal chances of taking 0 or 1, and Corr$\left( V_{1},V_{3}\right) =$Corr$\left( V_{1},V_{2}\right) =$Corr$\left( V_{2},V_{3}\right) =1/3$. $X_{j},$ $j=5,6,...,p_{n}$ are the same as those in DGP 9.\footnote{An example of $V_{1},V_{2},$ and $V_{3}$ is that $V_{j}=\mathbf{1(}\tilde{U}_{j}+\tilde{U}_{4}>1)$ for $j=1,2,$ and $3,$ where $\tilde{U}_{1},...,\tilde{U}_{4}$ are independent $U\left( 0,1\right) .$} Some simple calculation implies that $X_{4}$ is a hidden signal.
The key tuning parameter in this study is the threshold $\varsigma_{n}.$\ The CKP's Bonferroni correction strategy does not perform consistently well for our case. We do the following instead. We set \[ \varsigma_{n}=Cm_{n}\left[ \left( \log p_{n}\right) ^{1.1}+\left( \log m_{n}\right) ^{1.1}\right] \text{ and }\varsigma_{n}=C\left( \log p_{n}\right) ^{1.1} \] for continuous variables and binary variables, respectively. This satisfies Assumption (ref) when, in addition, Assumptions (ref) and (ref) hold (both $m_{n}$ and $p_{n}$ grow at the polynomial rate of $n$). In the small samples, we set this $\varsigma_{n}$ to have control over both $m_{n}$ and $p_{n}$. CKP suggest using a larger threshold, $\varsigma _{n}^{\ast}$, for subsequent stages, to improve the finite sample performance. We follow their lead and set a $\varsigma_{n}^{\ast}$ larger than $\varsigma_{n}$\ for subsequent stages. That is, we replace $\varsigma_{n}$ with $\varsigma_{n}^{\ast}$ for $\widehat{\mathcal{J}}_{l,\left( 2\right) },\widehat{\mathcal{J}}_{l,\left( 3\right) },...,\widehat{\mathcal{J} }_{l,\left( k\right) }$ in Section (ref). The reason is that the chance of including a noise variable increases quickly for subsequent stages, and we need a larger $\varsigma_{n}^{\ast}$ for the OCMT to conclude more easily. We set $\varsigma_{n}^{\ast}=3\varsigma_{n}$ for continuous variables. Note that CKP take the threshold for later stages to be twice the threshold for the first stage, and their test is based on $t$-statistics. Since ours is based on chi-squared statistics, equivalently, we should set $\varsigma_{n}^{\ast}$ to be $4\varsigma_{n}.$ However, some small-scale experiments suggest that setting $\varsigma_{n}^{\ast}=3\varsigma_{n}$ can yield better small sample performance, probably because our model is nonparametric and our test statistic varies more than the parametric counterpart. Note that this change does not affect our theoretical results because $\varsigma_{n}^{\ast}$ is proportional to $\varsigma_{n}$. For binary additive components, we continue to follow CKP and set $\varsigma_{n}^{\ast }=4\varsigma_{n}$ because they enter the model linearly.
Another important tuning parameter is $m_{n},$ which is critical for the sieve estimation. The optimal choice of $m_{n}$ has been studied extensively in the literature. Popular ways to choose the value of $m_{n}$ include cross validation, Akaike information criterion (AIC), and BIC. We refer the readers to Chen_handbook and Hansen2014 for a review on this important issue. For $m_{n},$ we simply set $m_{n}=\left\lfloor n^{1/4}\right\rfloor +1,$ where $\left\lfloor \cdot\right\rfloor $ is the floor operator$.$ This $m_{n}$ satisfies Assumption (ref), if $d>3/2$. Other choices of sieve terms such as $m_{n}=\left\lfloor n^{1/4}\right\rfloor +2$\ are also considered in the simulations, and they yield similar results.
For the $C$ in $\varsigma_{n}=Cm_{n}\left[ \left( \log p_{n}\right) ^{1.1}+\left( \log m_{n}\right) ^{1.1}\right] $ or $\varsigma_{n}=C\left( \log p_{n}\right) ^{1.1},$\ we test $C$ in the range of $0.5$ to $2.5$, specifically, $0.5,0.6,...,2.5$. We determine the value of $C$ by minimizing the following BIC:
where RSS$\left( C\right) $ denotes the residual sum of squares from the post-OCMT ordinary least squares regression of $Y$ on $P^{m_{n}}\left[ \boldsymbol{Z}_{\left( \hat{k}_{s}\right) }\right] $, where $\boldsymbol{Z} _{\left( \hat{k}_{s}\right) }$ are the variables selected by the OCMT with $C$ in use. The tuning parameters for the adaptive group Lasso\ are selected as in Section (ref). Another popular way to choose $C$ is cross validation. As seen from the results, the BIC works well for our procedure. In light of this, we will not pursue the procedure with cross validation.
For an easy reference, we present the implementation details in Appendix (ref).
For comparison, we consider the adaptive group Lasso by HuangEtal2010, which is designed for component selection in the nonparametric additive model. It is a two-step approach---the first step is the usual group Lasso and the second is the adaptive group Lasso with initial estimates from the first step. We select tuning parameters $\lambda_{n1}$ in step 1 and $\lambda_{n2}$ in step 2 ($\lambda_{n1}$ and $\lambda_{n2}$ are in the notations of HuangEtal2010) by BIC as well. Specifically, we set $\lambda _{n1}=\lambda_{j}$ with \[ \lambda_{j}=\exp\left\{ \log\left( \lambda_{\max}\right) +\left[ \log\left( \lambda_{\min}\right) -\log\left( \lambda_{\max}\right) \right] \frac{j}{30}\right\} \] for $j=0,1,...,30,$ where $\lambda_{\min}=\max\left\{ 0.05,10^{-5}\left\Vert \boldsymbol{y}\right\Vert \right\} $ and $\lambda_{\max}=0.5\left\Vert \boldsymbol{y}\right\Vert .$ We calculate BIC$_{j}$ based on the estimation using $\lambda_{n1}=\lambda_{j},$ and we select the $\hat{\lambda}_{n1}$ that minimizes the BIC$_{j}$ among $j=0,1,...,30.$ We set the estimates in the first step as the estimates using $\hat{\lambda}_{n1}.$ For the second step, we set $\lambda_{n2}=\lambda_{j}$ for $j=0,1,...,30$ with the initial estimates as the estimates from the first step using $\hat{\lambda}_{n1}$. We again calculate BIC$_{j}$ based on the estimation using $\lambda_{n2} =\lambda_{j},$ and we select the $\hat{\lambda}_{n2}$ that minimizes the BIC$_{j}$ among $j=0,1,...,30.$ The final estimates are the adaptive group Lasso estimates using $\hat{\lambda}_{n2}$. Note that the variables not selected in the first step are not included in the second step, because the penalty for those variables is infinity in the second step. We find the solutions in both steps through the block coordinate descent algorithm (i.e., the \textquotedblleft shooting\textquotedblright\ algorithm).\ For the algorithm details, see WuLange2008.\footnote{One implementation in MATLAB can be found at: https://publish.illinois.edu/xiaohuichen/code/group-lasso-shooting/.}
In addition, we compare our method with the random forest regression and bagging, both of which are commonly-used machine learning methods. Since there is no variable selection criterion in random forest, we focus on the comparison of out-sample forecasting performance.
We consider the combinations of $n=200$ or $400$ and $p_{n}=100$, $200$, or $1,000$ for each DGP. All results are based on 1,000 replications. We report the results for five different methods. The first and second methods are our post-OCMT procedure (Steps 1--6 in Appendix (ref), denoted as \textquotedblleft POST--OCMT\textquotedblright) and OCMT procedure (Steps 1-- 5 in Appendix (ref), denoted as \textquotedblleft OCMT\textquotedblright), respectively. The third method is the adaptive group Lasso (denoted as \textquotedblleft AGLASSO\textquotedblright), and the last two methods are bagging (denoted as \textquotedblleft BAGGING\textquotedblright) and random forest (denoted as \textquotedblleft RF\textquotedblright), respectively.
Following the literature, we report the mean number of variables selected (NV), the true positive rates (TPR), the false positive rates (FPR), the false discovery rates (FDR), the percentage of correct selection (CS) for the first three methods\footnote{That is, we precisely uncover the true model in ((ref)) using all the true signals, but not any pseudo-signals or noise variables.}, and the out-of-sample root mean squared forecast errors (RMSFE) for all the five methods.\footnote{The forecasts are based on the post-selection estimates from each method. In each replication, we additionally independently generate 200 observations. Then, we calculate the root mean square errors (RMSE) of the difference between the forecasts and $Y$ for the new observations. RMSFE is the average of those RMSE for 1,000 replications.} We report the average number of stages (STEP) for our OCMT procedure only.
To save space, we report the results in the main body of this paper for only DGPs 1 and 6 in Tables (ref) and (ref), respectively. These two designs correspond to the simplest and the most complicated designs in the simulation, respectively. Results for DGPs 2--5 and 7--10 are provided in Tables (ref) to (ref) in Appendix (ref). To showcase the advantage of the post-OCMT procedure compared to the OCMT procedure only, we also consider a linear DGP (DGP 11) in the online Appendix (ref) with results reported in Table (ref).
We can make several observations. First, POST--OCMT performs the best among the three methods in most cases. POST--OCMT outperforms OCMT due to its ability to eliminate pseudo-signals or noise variables.\ POST--OCMT even outperforms OCMT slightly for DGPs without any pseudo-signals, especially for $n=200$. This result confirms the necessity of conducting the additional post-OCMT step. Moreover, the additional computation time for POST--OCMT compared with OCMT can be almost ignored because the number of candidate variables after OCMT is very small. Second, AGLASSO performs almost the same as our procedures for DGP 2. Note that DGP 2 contains four independent signals, two pseudo-signals, and no hidden signals. Then, this result is not surprising, because Lasso performs very well for independent signals, and the biggest challenge for our procedure is the possible presence of pseudo-signals. When the signals are correlated in DGP 4, AGLASSO performs less well and is outperformed by POST--OCMT. Third, the performance of all methods improves as $n$ increases from $200$ to $400$. Fourth, OCMT is very successful at picking up hidden signals, especially for $n=400$; see, for example, the results for DGPs 3, 4, 6, and 8. Fifth, our methods perform well in the presence of heteroskedastic errors for DGP 7 and 8. Sixth, the CS of POST--OCMT performs well even for the complicated DGPs 6 and 8 at $n=400$, whereas AGLASSO performs poorly in terms of CS for these two DGPs. Seventh, the number of stages for OCMT basically confirms our theoretical findings. For example, for DGP 1, the mean number of stages is slightly more than 1 for $n=200$, and is 1 for $n=400$. In the presence of hidden signals, the mean number of stages is approximately 2 for $n=400$ for DGPs\ 3, 4, 6, 8, and 10. Note the mean number of stages is approximately 2 for $n=400$ for DGP\ 5 with correlated signals and no hidden signals. The reason is that the correlation makes the net effect of $X_{1}$ on $Y$ very small and $X_{1}$ almost behaves like a \textquotedblleft hidden signal\textquotedblright. This result also confirms the necessity of using multiple stages instead of a single-stage procedure. Eighth, our method also works well for DGPs 9 and 10 that mimic the application. An additional remark is that our procedure can be implemented fast and is much faster than AGLASSO. For example, when $p=1,000,$ our procedure took less than half a minute for one replication on average, whereas the AGLASSO took hours for one replication. Finally, the first three procedures have smaller RMSFE than Bagging and RF in almost all scenarios and thus have better out-of-sample forecasting performance.\footnote{Note that Bagging has slightly better out-of-sample forecasting performance than RF. This is reasonable given that our DGPs have a finite number of signal variables. In RF, many decision splits do not improve predictive accuracy because they rely solely on noise variables to generate the trees.} This is because they explicitly utilize the information (in the form of an additive function) underlying the DGPs while BAGGING and RF do not. In addition, POST--OCMT delivers the best out-of-sample forecasting performance among all procedures.
To summarize, our methods perform well in small samples, and we view it as a useful alternative to existing methods in the literature.
In this section, we apply our method to a dataset extracted from RUMiC.\footnote{RUMiC consists of three parts: the Urban Household Survey, the Rural Household Survey, and the Migrant Household Survey. A group of researchers at the Australian National University, the University of Queensland, and the Beijing Normal University initiated this survey. The Institute for the Study of Labor (IZA) supported it and provides the Scientific Use Files. RUMiC had financial support from the Australian Research Council, the Australian Agency for International Development, the Ford Foundation, IZA, and the Chinese Foundation of Social Sciences. More information on the survey can be found at https://datasets.iza.org/dataset/58/longitudinal-survey-on-rural-urban-migration-in-china.} The survey studied immigrants or workers moving from the rural areas of China to its big cities. The survey asked interviewees (immigrants or workers) a wide range of questions. For the detailed design of the survey and other information, including on the construction of each variable, see the survey website and Gongetal2008. Currently, the 2008 wave data are publicly available.
Economic reforms since the late 1970s have brought significant changes to China's economy. The government began relaxing its policy on population mobility in the early 1980s. Gradually, peasants were allowed to leave villages and work in big cities to earn higher incomes. Most migrant workers may leave their spouses, children, or parents behind in their hometowns\, who may need their financial support. This situation results in monetary transfers, that is, remittances, from migrant workers to their family. In the context of migration, family, and economic development, remittances are not only an income source for recipients, but also reflect intrafamilial relationships. Remittances clearly represent a dimension of family ties and demonstrate high degrees of interaction between migrants and families at home. In addition, remittances from the rural migrant workers also contribute significantly to China's agricultural productivity (c.f., Rozelle_at_al.1999). For these reasons, it has long been of interest to model remittances to families or relatives in the hometown; see Li2001 and Cai2003, among others.
In this application, we take remittance as the dependent variable ($Y$); we focus on the dataset from Guangdong Province (Guangzhou, Dongguan, and Shenzhen cities) in the 2008 survey wave, and keep 78 covariates\ from the dataset.\footnote{We keep covariates with relatively fewer missing observations.} After dropping observations with missing information, the number of observations is 456. We provide the definitions of the dependent variable and covariates, and the associated summary statistics, in Tables (ref) and (ref), respectively. We report the original labels of all covariates in the survey in the first column of Table (ref) for reference. Among the 78 covariates, there are some continuous variables with most observations as 0 (Panel C in Table (ref)), and some discrete variables with very limited support (Panel D in Table (ref)). For those variables, the design matrix of the sieves generated are either singular or close to singular. For this reason, we add those variables linearly into the model and treat them the same as\ dummy variables (Panel E in Table (ref)) for modeling. Consequently, the way we fit the dataset resembles the approach we used for DGPs 9 and 10 in the simulation. We explain the reason our theoretical results continue to hold in this situation in the simulation section (DGP 9). We take natural logarithms for the $Y$ and continuous $X$ variables to offset the effect of outliers; otherwise, the forecast can easily take some\ extreme values.
We randomly select 400 observations as the training sample, and the remaining 56 observations as the test sample. The number of sieve terms is set as $m_{n}=\left\lfloor 400^{1/4}\right\rfloor +1=5$ for the continuous variables in Panel B of Table (ref)). We set $\varsigma_{n}=Cm_{n}\left[ \left( \log p_{n}\right) ^{1.1}+\left( \log m_{n}\right) ^{1.1}\right] $ and $\varsigma_{n}^{\ast}=3\varsigma_{n}$ for continuous variables, and $\varsigma_{n}=C\left( \log p_{n}\right) ^{1.1}$ and $\varsigma_{n}^{\ast }=4\varsigma_{n}$ for terms entering the model linearly (see Section (ref) for the reason). We set $C$ in the range of $0.5$ to $2.5$, specifically, $0.5,0.6,...,2.5$, and choose $C$ to minimize the BIC for the model selection$.$ The competing methods are the group Lasso (labelled as GLASSO) and the adaptive group Lasso (AGLASSO). The tuning parameters for GLASSO and AGLASSO\ are selected as in Section (ref). We evaluate the performance of all methods based on the RMSFE of the test dataset using the fitted models from different methods. We independently repeat the above procedure 100 times.
We report the results in terms of the out-of-sample RMSFE in Tables (ref) and (ref), with AGLASSO as the benchmark. Table (ref) shows that OCMT and POST-OCMT outperform AGLASSO in the majority of cases. Of course, our methods do not outperform AGLASSO all the time. To highlight the necessity of the multiple stages, we also report the results of the one-stage procedure (we selected the tuning parameter also by minimizing the BIC) in Table (ref). The OCMT stops at the second stage for all 100 cases. We normalize the average RMSFE of AGLASSO as 1.\ Notably, the average RMSFE of our methods are lower than that of AGLASSO. The OCMT also improves the RMSFE over the one-stage procedure, possibly owing to some hidden signals uncovered by our multiple-stage procedure. We report the frequencies of variables (out of 100 cases) selected by all methods in Table (ref). It appears that G102 (monthly income) along with G133 (gifts to others, including parents) and G137 (education cost for left-behind children) contribute most to the model, as shown by all methods in general. We note that AGLASSO tends to select more variables than our methods, on average, which was also CKP's finding in their application. The \textquotedblleft over-fitting\textquotedblright\ is the main cause of the relative inferior performance of AGLASSO.
In this paper, we examine the one-covariate-at-a-time multiple testing approach to model selection in additive models. The properties of the TPR, FPR, and FDR of our approach are established based on some asymptotic probability bounds of Type-I and II errors. The simulation experiments and one application on the RUMiC dataset showcase excellent small-sample properties of our methods. Just as stated by CKP for linear models, we view our approach as a useful alternative to the model selection methods for additive models in the literature.