Variable selection in nonparametric additive models

Introduction

Let $(Y_{i},{\mathbf{X}}_{i}),i=1,\ldots,n$, be random vectors that are independently and identically distributed as $(Y,{\mathbf{X}})$, where $Y$ is a response variable and ${\mathbf{X}}=(X_{1},\ldots,X_{p})^{\prime}$ is a $p$-dimensional covariate vector. Consider the nonparametric additive model

where $\mu$ is an intercept term, $X_{ij}$ is the $j$th component of $X_{i}$, the $f_{j}$’s are unknown functions, and $\varepsilon_{i}$ is an unobserved random variable with mean zero and finite variance $\sigma^{2}$. Suppose that some of the additive components $f_{j}$ are zero. The problem addressed in this paper is to distinguish the nonzero components from the zero components and estimate the nonzero components. We allow the possibility that $p$ is larger than the sample size $n$, which we represent by letting $p$ increase as $n$ increases. We propose a penalized method for variable selection in (1) and show that the proposed method can correctly select the nonzero components with high probability.

There has been much work on penalized methods for variable selection and estimation with high-dimensional data. Methods that have been proposed include the bridge estimator [Frank and Friedman (1993), Huang, Horowitz and Ma (2008)]; least absolute shrinkage and selection operator or Lasso [Tibshirani (1996)], the smoothly clipped absolute deviation (SCAD) penalty [Fan and Li (2001), Fan and Peng (2004)], and the minimum concave penalty [Zhang (2010)]. Much progress has been made in understanding the statistical properties of these methods. In particular, many authors have studied the variable selection, estimation and prediction properties of the Lasso in high-dimensional settings. See, for example, Meinshausen and Bühlmann (2006), Zhao and Yu (2006), Zou (2006), Bunea, Tsybakov and Wegkamp (2007), Meinshausen and Yu (2009), Huang, Ma and Zhang (2008), van de Geer (2008) and Zhang and Huang (2008), among others. All these authors assume a linear or other parametric model. In many applications, however, there is little a priori justification for assuming that the effects of covariates take a linear form or belong to any other known, finite-dimensional parametric family. For example, in studies of economic development, the effects of covariates on the growth of gross domestic product can be nonlinear. Similarly, there is evidence of nonlinearity in the gene expression data used in the empirical example in Section 5.

There is a large body of literature on estimation in nonparametric additive models. For example, Stone (1985, 1986) showed that additive spline estimators achieve the same optimal rate of convergence for a general fixed $p$ as for $p=1$. Horowitz and Mammen (2004) and Horowitz, Klemelä and Mammen (2006) showed that if $p$ is fixed and mild regularity conditions hold, then oracle-efficient estimates of the $f_{j}$’s can be obtained by a two-step procedure. Here, oracle efficiency means that the estimator of each $f_{j}$ has the same asymptotic distribution that it would have if all the other $f_{j}$’s were known. However, these papers do not discuss variable selection in nonparametric additive models.

Antoniadis and Fan (2001) proposed a group SCAD approach for regularization in wavelets approximation. Zhang et al. (2004) and Lin and Zhang (2006) have investigated the use of penalization methods in smoothing spline ANOVA with a fixed number of covariates. Zhang et al. (2004) used a Lasso-type penalty but did not investigate model-selection consistency. Lin and Zhang (2006) proposed the component selection and smoothing operator (COSSO) method for model selection and estimation in multivariate nonparametric regression models. For fixed $p$, they showed that the COSSO estimator in the additive model converges at the rate $n^{-d/(2d+1)}$, where $d$ is the order of smoothness of the components. They also showed that, in the special case of a tensor product design, the COSSO correctly selects the nonzero additive components with high probability. Zhang and Lin (2006) considered the COSSO for nonparametric regression in exponential families.

Several other recent papers have also considered variable selection in nonparametric models. For example, Wang, Chen and Li (2007) and Wang and Xia (2008) considered the use of group Lasso and SCAD methods for model selection and estimation in varying coefficient models with a fixed number of coefficients and covariates. Bach (2007) applies what amounts to the group Lasso to a nonparametric additive model with a fixed number of covariates. He established model selection consistency under conditions that are considerably more complicated than the ones we require for a possibly diverging number of covariates.

In this paper, we propose to use the adaptive group Lasso for variable selection in (1) based on a spline approximation to the nonparametric components. With this approximation, each nonparametric component is represented by a linear combination of spline basis functions. Consequently, the problem of component selection becomes that of selecting the groups of coefficients in the linear combinations. It is natural to apply the group Lasso method, since it is desirable to take into the grouping structure in the approximating model. To achieve model selection consistency, we apply the group Lasso iteratively as follows. First, we use the group Lasso to obtain an initial estimator and reduce the dimension of the problem. Then we use the adaptive group Lasso to select the final set of nonparametric components. The adaptive group Lasso is a simple generalization of the adaptive Lasso [Zou (2006)] to the method of the group Lasso [Yuan and Lin (2006)]. However, here we apply this approach to nonparametric additive modeling.

We assume that the number of nonzero $f_{j}$’s is fixed. This enables us to achieve model selection consistency under simple assumptions that are easy to interpret. We do not have to impose compatibility or irrepresentable conditions, nor do we need to assume conditions on the eigenvalues of certain matrices formed from the spline basis functions. We show that the group Lasso selects a model whose number of components is bounded with probability approaching one by a constant that is independent of the sample size. Then using the group Lasso result as the initial estimator, the adaptive group Lasso selects the correct model with probability approaching 1 and achieves the optimal rate of convergence for nonparametric estimation of an additive model.

The remainder of the paper is organized as follows. Section 2 describes the group Lasso and the adaptive group Lasso for variable selection in nonparametric additive models. Section 3 presents the asymptotic properties of these methods in “large $p$, small $n$” settings. Section 4 presents the results of simulation studies to evaluate the finite-sample performance of these methods. Section 5 provides an illustrative application, and Section 6 includes concluding remarks. Proofs of the results stated in Section 3 are given in the Appendix.

Adaptive group Lasso in nonparametric additive models

We describe a two-step approach that uses the group Lasso for variable selection based on a spline representation of each component in additive models. In the first step, we use the standard group Lasso to achieve an initial reduction of the dimension in the model and obtain an initial estimator of the nonparametric components. In the second step, we use the adaptive group Lasso to achieve consistent selection.

There exists a normalized B-spline basis $\{\phi_{k},1\leq k\leq m_{n}\}$ for $\mathcal{S}_{n}$, where $m_{n}\equiv K_{n}+l$ [Schumaker (1981)]. Thus, for any $f_{nj}\in\mathcal{S}_{n}$, we can write

Under suitable smoothness assumptions, the $f_{j}$’s can be well approximated by functions in $\mathcal{S}_{n}$. Accordingly, the variable selection method described in this paper is based on the representation (2).

where $\lambda_{n}$ is a penalty parameter. We study the estimators that minimize $L_{n}(\mu,\bolds\beta_{n})$ subject to the constraints

For simplicity and without causing confusion, we simply write $\psi_{k}(x)=\psi_{jk}(x)$. Define

So, $Z_{ij}$ consists of values of the (centered) basis functions at the $i$th observation of the $j$th covariate. Let ${\mathbf{Z}}_{j}=(Z_{1j},\ldots,Z_{nj})^{\prime}$ be the $n\times m_{n}$ “design” matrix corresponding to the $j$th covariate. The total “design” matrix is ${\mathbf{Z}}=({\mathbf{Z}}_{1},\ldots,{\mathbf{Z}}_{p})$. Let ${\mathbf{Y}}=(Y_{1}-\overline{Y},\ldots,Y_{n}-\overline{Y})^{\prime}$. With this notation, we can write

Here, we have dropped $\mu$ in the argument of $L_{n}$. With the centering, $\widehat{\mu}=\overline{Y}$. Then minimizing (3) subject to (4) is equivalent to minimizing (6) with respect to $\bolds\beta_{n}$, but the centering constraints are not needed for (6).

We now describe the two-step approach to component selection in the nonparametric additive model (1).

Step 1. Compute the group Lasso estimator. Let

This objective function is the special case of (6) that is obtained by setting $w_{nj}=1$, $1\leq j\leq p$. The group Lasso estimator is $\widetilde{\bolds\beta}_{n}\equiv\widetilde{\bolds\beta}_{n}(\lambda_{n1})=\arg\min_{\bolds\beta_{n}}L_{n1}(\bolds\beta_{n};\lambda_{n1}).$

Step 2. Use the group Lasso estimator $\widetilde{\bolds\beta}_{n}$ to obtain the weights by setting

The adaptive group Lasso objective function is

Here, we define $0\cdot\infty=0$. Thus, the components not selected by the group Lasso are not included in Step 2. The adaptive group Lasso estimator is $\widehat{\bolds\beta}_{n}\equiv\widehat{\bolds\beta}_{n}(\lambda_{n2})=\arg\min_{\bolds\beta_{n}}L_{n2}(\bolds\beta_{n};\lambda_{n2})$. Finally, the adaptive group Lasso estimators of $\mu$ and $f_{j}$ are

Main results

This section presents our results on the asymptotic properties of the estimators defined in Steps 1 and 2 of Section 2.

Let $k$ be a nonnegative integer, and let $\alpha\in(0,1]$ be such that $d=k+\alpha>0.5$. Let $\mathcal{F}$ be the class of functions $f$ on $$whose$k$th derivative$f^{(k)}$exists and satisfies a Lipschitz condition of order$\alpha$:

In (1), without loss of generality, suppose that the first $q$ components are nonzero, that is, $f_{j}(x)\neq 0,1\leq j\leq q$, but $f_{j}(x)\equiv 0,q+1\leq j\leq p$. Let $A_{1}=\{1,\ldots,q\}$ and $A_{0}=\{q+1,\ldots,p\}$. Define $\|f\|_{2}=[\int_{a}^{b}f^{2}(x)\,dx]^{1/2}$ for any function $f$, whenever the integral exists.

(A1) The number of nonzero components $q$ is fixed and there is a constant $c_{f}>0$ such that ${\min_{1\leq j\leq q}}\|f_{j}\|_{2}\geq c_{f}$.

(A4) The covariate vector $X$ has a continuous density and there exist constants $C_{1}$ and $C_{2}$ such that the density function $g_{j}$ of $X_{j}$ satisfies $0<C_{1}\leq g_{j}(x)\leq C_{2}<\infty$ on $[a,b]$ for every $1\leq j\leq p$.

In this section, we consider the selection and estimation properties of the group Lasso estimator. Define ${\widetilde{A}}_{1}=\{j\dvtx\|\widetilde{\bolds\beta}_{nj}\|_{2}\neq 0,1\leq j\leq p\}$. Let $|A|$ denote the cardinality of any set $A\subseteq\{1,\ldots,p\}$.

Suppose that (A1) to (A4) hold and $\lambda_{n1}\geq C\sqrt{n\log(pm_{n})}$ for a sufficiently large constant $C$.

With probability converging to 1, $|{\widetilde{A}}_{1}|\leq M_{1}|A_{1}|=M_{1}q$ for a finite constant $M_{1}>1$.

If ${m_{n}^{2}\log(pm_{n})}/{n}\rightarrow 0$ and $(\lambda_{n1}^{2}m_{n})/n^{2}\rightarrow 0$ as $n\rightarrow\infty$, then all the nonzero $\bolds\beta_{nj},1\leq j\leq q$, are selected with probability converging to one.

Part (i) of Theorem 1 says that, with probability approaching 1, the group Lasso selects a model whose dimension is a constant multiple of the number of nonzero additive components $f_{j}$, regardless of the number of additive components that are zero. Part (ii) implies that every nonzero coefficient will be selected with high probability. Part (iii) shows that the difference between the coefficients in the spline representation of the nonparametric functions in (1) and their estimators converges to zero in probability. The rate of convergence is determined by four terms: the stochastic error in estimating the nonparametric components (the first term) and the intercept $\mu$ (the second term), the spline approximation error (the third term) and the bias due to penalization (the fourth term).

Let $\widetilde{f}_{nj}(x)=\sum_{j=1}^{m_{n}}\widetilde{\beta}_{jk}\psi(x),1\leq j\leq p$. The following theorem is a consequence of Theorem 1.

Suppose that (A1) to (A4) hold and that $\lambda_{n1}\geq\break C\sqrt{n\log(pm_{n})}$ for a sufficiently large constant $C$. Then:

Let ${\widetilde{A}}_{f}=\{j\dvtx\|\widetilde{f}_{nj}\|_{2}>0,1\leq j\leq p\}$. There is a constant $M_{1}>1$ such that, with probability converging to 1, $|{\widetilde{A}}_{f}|\leq M_{1}q$.

If $(m_{n}\log(pm_{n}))/n\rightarrow 0$ and $(\lambda_{n1}^{2}m_{n})/n^{2}\rightarrow 0$ as $n\rightarrow\infty$, then all the nonzero additive components $f_{j},1\leq j\leq q$, are selected with probability converging to one.

where ${\widetilde{A}}_{2}=A_{1}\cup{\widetilde{A}}_{1}$.

Thus, under the conditions of Theorem 2, the group Lasso selects all the nonzero additive components with high probability. Part (iii) of the theorem gives the rate of convergence of the group Lasso estimator of the nonparametric components.

For any two sequences $\{a_{n},b_{n},n=1,2,\ldots\}$, we write $a_{n}\asymp b_{n}$ if there are constants $0<c_{1}<c_{2}<\infty$ such that $c_{1}\leq a_{n}/b_{n}\leq c_{2}$ for all $n$ sufficiently large.

We now state a useful corollary of Theorem 2.

Suppose that (A1) to (A4) hold. If $\lambda_{n1}\asymp\sqrt{n\log(pm_{n})}$ and $m_{n}\asymp n^{1/(2d+1)}$, then:

If $n^{-2d/(2d+1)}\log(p)\rightarrow 0$ as $n\rightarrow\infty$, then with probability converging to one, all the nonzero components $f_{j},1\leq j\leq q$, are selected and the number of selected components is no more than $M_{1}q$.

For the $\lambda_{n1}$ and $m_{n}$ given in Corollary 1, the number of zero components can be as large as $\exp(o(n^{2d/(2d+1)}))$. For example, if each $f_{j}$ has continuous second derivative ($d=2$), then it is $\exp(o(n^{4/5}))$, which can be much larger than $n$.

2 Selection consistency of the adaptive group Lasso

We now consider the properties of the adaptive group Lasso. We first state a general result concerning the selection consistency of the adaptive group Lasso, assuming an initial consistent estimator is available. We then apply to the case when the group Lasso is used as the initial estimator. We make the following assumptions.

(B1) The initial estimators $\widetilde{\bolds\beta}_{nj}$ are $r_{n}$-consistent at zero:

and there exists a constant $c_{b}>0$ such that

where $b_{n1}={\min_{j\in A_{1}}}\|\bolds\beta_{nj}\|_{2}$.

(B2) Let $q$ be the number of nonzero components and $s_{n}=p-q$ be the number of zero components. Suppose that:

We state condition (B1) for a general initial estimator, to highlight the point that the availability of an $r_{n}$-consistent estimator at zero is crucial for the adaptive group Lasso to be selection consistent. In other words, any initial estimator satisfying (B1) will ensure that the adaptive group Lasso (based on this initial estimator) is selection consistent, provided that certain regularity conditions are satisfied. We note that it follows immediately from Theorem 1 that the group Lasso estimator satisfies (B1). We will come back to this point below.

For $\widehat{\bolds\beta}_{n}\equiv(\widehat{\bolds\beta}_{n1}^{\prime},\ldots,\widehat{\bolds\beta}_{np}^{\prime})^{\prime}$ and $\bolds\beta_{n}\equiv(\bolds\beta_{n1}^{\prime},\ldots,\bolds\beta_{np}^{\prime})^{\prime}$, we say $\widehat{\bolds\beta}_{n}=_{0}\bolds\beta_{n}$ if $\operatorname{sgn}_{0}(\|\widehat{\bolds\beta}_{nj}\|)=\operatorname{sgn}_{0}(\|\bolds\beta_{nj}\|),1\leq j\leq p$, where $\operatorname{sgn}_{0}(|x|)=1$ if $|x|>0$ and $=0$ if $|x|=0$.

Suppose that conditions (B1), (B2) and (A1)–(A4) hold. Then:

This theorem is concerned with the selection and estimation properties of the adaptive group Lasso in terms of $\widehat{\bolds\beta}_{n}$. The following theorem states the results in terms of the estimators of the nonparametric components.

Suppose that conditions (B1), (B2) and (A1)–(A4) hold. Then:

Part (i) of this theorem states that the adaptive group Lasso can consistently distinguish nonzero components from zero components. Part (ii) gives an upper bound on the rate of convergence of the estimator.

We now apply the above results to our proposed procedure described in Section 2, in which we first obtain the the group Lasso estimator and then use it as the initial estimator in the adaptive group Lasso.

By Theorem 1, if $\lambda_{n1}\asymp\sqrt{n\log(pm_{n})}$ and $m_{n}\asymp n^{1/(2d+1)}$ for $d\geq 1$, then the group Lasso estimator satisfies (B1) with $r_{n}\asymp n^{d/(2d+1)}/\sqrt{\log(pm_{n})}$. In this case, (B2) simplifies to

We summarize the above discussion in the following corollary.

Let the group Lasso estimator $\widetilde{\bolds\beta}_{n}\equiv\widetilde{\bolds\beta}_{n}(\lambda_{n1})$ with $\lambda_{n1}\asymp\sqrt{n\log(pm_{n})}$ and $m_{n}\asymp n^{1/(2d+1)}$ be the initial estimator in the adaptive group Lasso. Suppose that the conditions of Theorem 1 hold. If $\lambda_{n2}\leq O(n^{1/2})$ and satisfies (7), then the adaptive group Lasso consistently selects the nonzero components in (1), that is, part (i) of Theorem 4 holds. In addition,

This corollary follows directly from Theorems 1 and 4. The largest $\lambda_{n2}$ allowed is $\lambda_{n2}=O(n^{1/2})$. With this $\lambda_{n2}$, the first equation in (6) is satisfied. Substitute it into the second equation in (6), we obtain $p=\exp(o(n^{2d/(2d+1)}))$, which is the largest $p$ permitted and can be larger than $n$. Thus, under the conditions of this corollary, our proposed adaptive group Lasso estimator using the group Lasso as the initial estimator is selection consistent and achieves optimal rate of convergence even when $p$ is larger than $n$. Following model selection, oracle-efficient, asymptotically normal estimators of the nonzero components can be obtained by using existing methods.

Simulation studies

We use simulation to evaluate the performance of the adaptive group Lasso with regard to variable selection. The generating model is

Since $p$ can be larger than $n$, we consider two ways to select the penalty parameter, the BIC [Schwarz (1978)] and the EBIC [Chen and Chen (2008, 2009)]. The BIC is defined as

where $0\leq\nu\leq 1$ is a constant. We use $\nu=0.5$.

We have also considered two other possible ways of defining df: (a) using the trace of a linear smoother based on a quadratic approximation; (b) using the number of estimated nonzero components. We have decided to use the definition given above based on the results from our simulations. We note that the df for the group Lasso of Yuan and Lin (2006) requires an initial (least squares) estimator, which is not available when $p>n$. Thus, their df is not applicable to our problem.

In our simulation example, we compare the adaptive group Lasso with the group Lasso and ordinary Lasso. Here, the ordinary Lasso estimator is defined as the value that minimizes

This simple application of the Lasso does not take into account the grouping structure in the spline expansions of the components. The group Lasso and the adaptive group Lasso estimates are computed using the algorithm proposed by Yuan and Lin (2006). The ordinary Lasso estimates are computed using the Lars algorithms [Efron et al. (2004)]. The group Lasso is used as the initial estimate for the adaptive group Lasso.

We also compare the results from the nonparametric additive modeling with those from the standard linear regression model with Lasso. We note that this is not a fair comparison because the generating model is highly nonlinear. Our purpose is to illustrate that it is necessary to use nonparametric models when the underlying model deviates substantially from linear models in the context of variable selection with high-dimensional data and that model misspecification can lead to bad selection results.

where $f_{1}(t)=5t,f_{2}(t)=3(2t-1)^{2},f_{3}(t)=4{\sin(2\pi t)}/{(2-\sin(2\pi t))},$ $f_{4}(t)=6(0.1\sin(2\pi t)+0.2\cos(2\pi t)+0.3\sin(2\pi t)^{2}+0.4\cos(2\pi t)^{3}+0.5\sin(2\pi t)^{3}),$ and $f_{5}(t)=\cdots=f_{p}(t)=0$. Thus, the number of nonzero functions is $q=4$. This generating model is the same as Example 1 of Lin and Zhang (2006). However, here we use this model in high-dimensional settings. We consider the cases where $p=1000$ and three different sample sizes: $n=50,100$ and $200$. We use the cubic B-spline with six evenly distributed knots for all the functions $f_{k}$. The number of replications in all the simulations is 400.

The covariates are simulated as follows. First, we generate $w_{i1},\ldots,w_{ip},u_{i}$, $u^{\prime}_{i},v_{i}$ independently from $N(0,1)$ truncated to the interval $$,$i=1,\ldots,n$. Then we set$x_{ik}=(w_{ik}+tu_{i})/(1+t)$for$k=1,\ldots,4$and$x_{ik}=(w_{ik}+tv_{i})/(1+t)$for$k=5,\ldots,p$, where the parameter$t$controls the amount of correlation among predictors. We have$\operatorname{Corr}(x_{ik},x_{ij})=t^{2}/(1+t^{2})$,$1\leq j\leq 4$,$1\leq k\leq 4$, and$\operatorname{Corr}(x_{ik},x_{ij})=t^{2}/(1+t^{2})$,$4\leq j\leq p$,$4\leq k\leq p$, but the covariates of the nonzero components and zero components are independent. We consider$t=0,1$in our simulation. The signal to noise ratio is defined to be$sd(f)/sd(\epsilon)$. The error term is chosen to be$\epsilon_{i}\sim N(0,1.27^{2})$to give a signal-to-noise ratio (SNR)$3.11:1$. This value is the same as the estimated SNR in the real data example below, which is the square root of the ratio of the sum of estimated components squared divided by the sum of residual squared.

The results of 400 Monte Carlo replications are summarized in Table 1. The columns are the mean number of variables selected (NV), model error (ER), the percentage of replications in which all the correct additive components are included in the selected model (IN), and the percentage of replications in which precisely the correct components are selected (CS). The corresponding standard errors are in parentheses. The model error is computed as the average of $n^{-1}\sum_{i=1}^{n}[\hat{f}(x_{i})-f(x_{i})]^{2}$ over the 400 Monte Carlo replications, where $f$ is the true conditional mean function.

Table 1 shows that the adaptive group Lasso selects all the nonzero components (IN) and selects exactly the correct model (CS) more frequently than the other methods do. For example, with the BIC and $n=200$, the percentage of correct selections (CS) by the adaptive group Lasso ranges from 65.25% to 81%, which is much higher than the ranges 30–57.75% for the group Lasso and 12–15.75% for the ordinary Lasso. The adaptive group Lasso and group Lasso perform better than the ordinary Lasso in all of the experiments, which illustrates the importance of taking account of the group structure of the coefficients of the spline expansion. Correlation among covariates increases the difficulty of component selection, so it is not surprising that all methods perform better with independent covariates than with correlated ones. The percentage of correct selections increases as the sample size increases. The linear model with Lasso never selects the correct model. This illustrates the poor results that can be produced by a linear model when the true conditional mean function is nonlinear.

Table 1 also shows that the model error (ME) of the group Lasso is only slightly larger than that of the adaptive group Lasso. The models selected by the group Lasso nest and, therefore, have more estimated coefficients than the models selected by the adaptive group Lasso. Therefore, the group Lasso estimators of the conditional mean function have a larger variance and larger ME. The differences between the MEs of the two methods are small, however, because as can be seen from the NV column, the models selected by the group Lasso in our experiments have only slightly more estimated coefficients than the models selected by the adaptive group Lasso.

We now compare the adaptive group Lasso with the COSSO [Lin and Zhang (2006)]. This comparison is suggested to us by the Associate Editor. Because the COSSO algorithm only works for the case when $p$ is smaller than $n$, we use the same set-up as in Example 1 of Lin and Zhang (2006). In this example, the generating model is as in (8) with 4 nonzero components. Let $X_{j}=(W_{j}+tU)/(1+t)$, $j=1,\ldots,p$, where $W_{1},\ldots,W_{p}$ and $U$ are i.i.d. from $N(0,1)$, truncated to the interval $$. Therefore, corr$(X_{j},X_{k})=t^{2}/(1+t^{2})$for$j\neq k$. The random error term$\epsilon\sim N(0,1.32^{2})$. The SNR is 3:1. We consider three different sample sizes$n=50,100$or$200$and three different number of predictors$p=10,20$or$50$. The COSSO estimator is computed using the Matlab software which is publicly available at http://www4.stat.ncsu.edu/~hzhang/cosso.html.

The COSSO procedure uses either generalized cross-validation or 5-fold cross-validation. Based the simulation results of Lin and Zhang (2006) and our own simulations, the COSSO with 5-fold cross-validation has better selection performance. Thus, we compare the adaptive group Lasso with BIC or EBIC with the COSSO with 5-fold cross-validation. The results are given in Table 2. For independent predictors, when $n=200$ and $p=10,20$ or $50$, the adaptive group Lasso and COSSO have similar performance in terms of selection accuracy and model error. However, for smaller $n$ and larger $p$, the adaptive group Lasso does significantly better. For example, for $n=100$ and $p=50$, the percentage of correct selection for the adaptive group Lasso is 81–83%, but it is only 11% for the COSSO. The model error of the adaptive group Lasso is similar to or smaller than that of the COSSO. In several experiments, the model error of the COSSO is 2 to more than 7 times larger than that of the adaptive group Lasso. It is interesting to note that when $n=50$ and $p=20$ or $50$, the adaptive group Lasso still does a descent job in selecting the correct model, but the COSSO does poorly in these two cases. In particular, for $n=50$ and $p=50$, the COSSO did not select the exact correct model in all the simulation runs. For dependent predictors, the comparison is even mode favorable to the adaptive group Lasso, which performs significantly better than COSSO in terms of both model error and selection accuracy in all the cases.

Data example

We use the data set reported in Scheetz et al. (2006) to illustrate the application of the proposed method in high-dimensional settings. For this data set, 120 twelve-week old male rats were selected for tissue harvesting from the eyes and for microarray analysis. The microarrays used to analyze the RNA from the eyes of these animals contain over 31,042 different probe sets (Affymetric GeneChip Rat Genome 230 2.0 Array). The intensity values were normalized using the robust multi-chip averaging method [Irizzary et al. (2003)] method to obtain summary expression values for each probe set. Gene expression levels were analyzed on a logarithmic scale.

We are interested in finding the genes that are related to the gene TRIM32. This gene was recently found to cause Bardet–Biedl syndrome [Chiang et al. (2006)], which is a genetically heterogeneous disease of multiple organ systems including the retina. Although over 30,000 probe sets are represented on the Rat Genome 230 2.0 Array, many of them are not expressed in the eye tissue and initial screening using correlation shows that most probe sets have very low correlation with TRIM32. In addition, we are expecting only a small number of genes to be related to TRIM32. Therefore, we use 500 probe sets that are expressed in the eye and have highest marginal correlation in the analysis. Thus, the sample size is $n=120$ (i.e., there are 120 arrays from 120 rats) and $p=500$. It is expected that only a few genes are related to TRIM32. Therefore, this is a sparse, high-dimensional regression problem.

We use the nonparametric additive model to model the relation between the expression of TRIM32 and those of the 500 genes. We estimate model (1) using the ordinary Lasso, group Lasso, and adaptive group Lasso for the nonparametric additive model. To compare the results of the nonparametric additive model with that of the linear regression model, we also analyzed the data using the linear regression model with Lasso. We scale the covariates so that their values are between 0 and 1 and use cubic splines with six evenly distributed knots to estimate the additive components. The penalty parameters in all the methods are chosen using the BIC or EBIC as in the simulation study. Table 3 lists the probes selected by the group Lasso and the adaptive group Lasso, indicated by the check signs. Table 4 shows the number of variables, the residual sums of squares obtained with each estimation method. For the ordinary Lasso with the spline expansion, a variable is considered to be selected if any of the estimated coefficients of the spline approximation to its additive component are nonzero. Depending on whether BIC or EBIC is used, the group Lasso selects 16–17 variables, the adaptive group Lasso selects 15 variables and the ordinary Lasso with the spline expansion selects 94–97 variables, the linear model selects 8–14 variables. Table 4 shows that the adaptive group Lasso does better than the other methods in terms of residual sum of squares (RSS). We have also examined the plots (not shown) of the estimated additive components obtained with the group Lasso and the adaptive group Lasso, respectively. Most are highly nonlinear, confirming the need for taking into account nonlinearity.

In order to evaluate the performance of the methods, we use cross-validation and compare the prediction mean square errors (PEs). We randomly partition the data into 6 subsets, each set consisting of 20 observations. We then fit the model with 5 subsets as training set and calculate the PE for the remaining set which we consider as test set. We repeat this process 6 times, considering one of the 6 subsets as test set every time. We compute the average of the numbers of probes selected and the prediction errors of these 6 calculations. Then we replicate this process 400 times (this is suggested to us by the Associate Editor). Table 5 gives the average values over 400 replications. The adaptive group Lasso has smaller average prediction error than the group Lasso, the ordinary Lasso and the linear regression with Lasso. The ordinary Lasso selects far more probe sets than the other approaches, but this does not lead to better prediction performance. Therefore, in this example, the adaptive group Lasso provides the investigator a more targeted list of probe sets, which can serve as a starting point for further study.

It is of interest to compare the selection results from the adaptive group Lasso and the linear regression model with Lasso. The adaptive group Lasso and the linear model with Lasso select different sets of genes. When the penalty parameter is chosen with the BIC, the adaptive group Lasso selects 5 genes that are not selected by the linear model with Lasso. In addition, the linear model with Lasso selects 5 genes that are not selected by the adaptive group Lasso. When the penalty parameter is selected with the EBIC, the adaptive group Lasso selects 10 genes that are not selected by the linear model with Lasso. The estimated effects of many of the genes are nonlinear, and the Monte Carlo results of Section 4 show that the performance of the linear model with Lasso can be very poor in the presence of nonlinearity. Therefore, we interpret the differences between the gene selections of the adaptive group Lasso and the linear model with Lasso as evidence that the selections produced by the linear model are misleading.

Concluding remarks

In this paper, we propose to use the adaptive group Lasso for variable selection in nonparametric additive models in sparse, high-dimensional settings. A key requirement for the adaptive group Lasso to be selection consistent is that the initial estimator is estimation consistent and selects all the important components with high probability. In low-dimensional settings, finding an initial consistent estimator is relatively easy and can be achieved by many well-established approaches such as the additive spline estimators. However, in high-dimensional settings, finding an initial consistent estimator is difficult. Under the conditions stated in Theorem 1, the group Lasso is shown to be consistent and selects all the important components. Thus the group Lasso can be used as the initial estimator in the adaptive Lasso to achieve selection consistency. Following model selection, oracle-efficient, asymptotically normal estimators of the nonzero components can be obtained by using existing methods. Our simulation results indicate that our procedure works well for variable selection in the models considered. Therefore, the adaptive group Lasso is a useful approach for variable selection and estimation in sparse, high-dimensional nonparametric additive models.

Our theoretical results are concerned with a fixed sequence of penalty parameters, which are not applicable to the case where the penalty parameters are selected based on data driven procedures such as the BIC. This is an important and challenging problem that deserves further investigation, but is beyond the scope of this paper. We have only considered linear nonparametric additive models. The adaptive group Lasso can be applied to generalized nonparametric additive models, such as the generalized logistic nonparametric additive model and other nonparametric models with high-dimensional data. However, more work is needed to understand the properties of this approach in those more complicated models.

Appendix: Proofs

We first prove the following lemmas. Denote the centered versions of $\mathcal{S}_{n}$ by

where $\psi_{k}$’s are the centered spline bases defined in (5).

In particular, if we choose $m_{n}=O(n^{1/(2d+1)})$, then

By (A4), for $f\in\mathcal{F}$, there is an $f_{n}^{*}\in\mathcal{S}_{n}$ such that $\|f-f_{n}^{*}\|_{2}=O(m_{n}^{-d})$. Let $f_{n}=f_{n}^{*}-n^{-1}\sum_{i=1}^{n}f_{n}^{*}(X_{ij})$. Then $f_{n}\in\mathcal{S}^{0}_{nj}$ and $|f_{n}-f|\leq|f_{n}^{*}-f|+|P_{n}f_{n}^{*}|$, where $P_{n}$ is the empirical measure of i.i.d. random variables $X_{1j},\ldots,X_{nj}$. Consider

Here, we use the linear functional notation, for example, $Pf=\int fdP$, where $P$ is the probability measure of $X_{1j}$. For any $\varepsilon>0$, the bracketing number $N_{[\cdot]}(\varepsilon,\mathcal{S}_{nj}^{0},L_{2}(P))$ of $\mathcal{S}_{nj}^{0}$ satisfies $\log N_{[\cdot]}(\varepsilon,\mathcal{S}_{nj}^{0},L_{2}(P))\leq c_{1}m_{n}\log(1/\varepsilon)$ for some constant $c_{1}>0$ [Shen and Wong (1994), page 597]. Thus, by the maximal inequality; see, for example, van der Vaart (1998, page 288), $(P_{n}-P)f_{n}^{*}=O_{p}(n^{-1/2}m_{n}^{1/2})$. By (A4), $|P(f_{n}^{*}-f)|\leq C_{2}\|f_{n}^{*}-f\|_{2}=O(m_{n}^{-d})$ for some constant $C_{2}>0$. The lemma follows from the triangle inequality.

Suppose that conditions (A2) and (A4) hold. Let

and $T_{n}={\max_{1\leq j\leq p,1\leq k\leq m_{n}}}|T_{jk}|$. Then

where $C_{1}$ and $C_{2}$ are two positive constants. In particular, when $m_{n}\log(pm_{n})/\break n\rightarrow 0$,

Let $s_{njk}^{2}=\sum_{i=1}^{n}\psi_{k}^{2}(X_{ij})$. Conditional on $X_{ij}$’s, $T_{jk}$’s are sub-Gaussian. Let $s_{n}^{2}=\max_{1\leq j\leq p,1\leq k\leq m_{n}}s_{njk}^{2}$. By (A2) and the maximal inequality for sub-Gaussian random variables [van der Vaart and Wellner (1996), Lemmas 2.2.1 and 2.2.2],

where $C_{1}>0$ is a constant. By (A4) and the properties of B-splines,

for a constant $C_{2}>0$, for every $1\leq j\leq p$ and $1\leq k\leq m_{n}$. By (10),

By Lemma A.1 of van de Geer (2008), (10) and (11) imply

Therefore, by (12) and the triangle inequality,

Here, $\bolds\beta_{A}$ is an $|A|m_{n}\times 1$ vector and ${\mathbf{Z}}_{A}$ is an $n\times|A|m_{n}$ matrix. Let ${\mathbf{C}}_{A}={\mathbf{Z}}_{A}^{\prime}{\mathbf{Z}}_{A}/n.$ When $A=\{1,\ldots,p\}$, we simply write ${\mathbf{C}}={\mathbf{Z}}^{\prime}{\mathbf{Z}}/n$. Let $\rho_{\min}({\mathbf{C}}_{A})$ and $\rho_{\max}({\mathbf{C}}_{A})$ be the minimum and maximum eigenvalues of ${\mathbf{C}}_{A}$, respectively.

Let $m_{n}=O(n^{\gamma})$ where $0<\gamma<0.5$. Suppose that $|A|$ is bounded by a fixed constant independent of $n$ and $p$. Let $h\equiv h_{n}\asymp m_{n}^{-1}$. Then under (A3) and (A4), with probability converging to one,

where $c_{1}$ and $c_{2}$ are two positive constants.

Without loss of generality, suppose $A=\{1,\ldots,k\}$. Then ${\mathbf{Z}}_{A}=({\mathbf{Z}}_{1}$, $\ldots,{\mathbf{Z}}_{q})$. Let ${\mathbf{b}}=({\mathbf{b}}_{1}^{\prime},\ldots,{\mathbf{b}}_{q}^{\prime})^{\prime}$, where ${\mathbf{b}}_{j}\in R^{m_{n}}$. By Lemma 3 of Stone (1985),

for a certain constant $c_{3}>0$. By the triangle inequality,

Since ${\mathbf{Z}}_{A}{\mathbf{b}}={\mathbf{Z}}_{1}{\mathbf{b}}_{1}+\cdots+{\mathbf{Z}}_{q}{\mathbf{b}}_{q}$, the above two inequalities imply that

Let ${\mathbf{C}}_{j}=n^{-1}{\mathbf{Z}}_{j}^{\prime}{\mathbf{Z}}_{j}$. By Lemma 6.2 of Zhou, Shen and Wolf (1998),

Since ${\mathbf{C}}_{A}=n^{-1}{\mathbf{Z}}_{A}^{\prime}{\mathbf{Z}}_{A}$, it follows from (Appendix: Proofs) that

The lemma follows. {pf*}Proof of Theorem 1 The proof of parts (i) and (ii) essentially follows the proof of Theorem 2.1 of Wei and Huang (2008). The only change that must be made here is that we need to consider the approximation error of the regression functions by splines. Specifically, let $\bolds\xi_{n}=\bolds\varepsilon_{n}+\bolds\delta_{n}$, where $\bolds\delta_{n}=(\delta_{n1},\ldots,\delta_{nn})^{\prime}$ with $\delta_{ni}=\sum_{j=1}^{q_{n}}(f_{0j}(X_{ij})-f_{nj}(X_{ij}))$. Since $\|f_{0j}-f_{nj}\|_{2}=O(m_{n}^{-d})=O(n^{-d/(2d+1)})$ for $m_{n}=n^{1/(2d+1)}$, we have

for some constant $C_{1}>0$. For any integer $t$, let

where $V_{A}(S_{A})=\bolds\xi_{n}^{\prime}({\mathbf{Z}}_{A}({\mathbf{Z}}_{A}^{\prime}{\mathbf{Z}}_{A})^{-1}\bar{S}_{A}-(I-P_{A})X\bolds\beta$ for $N(A)=q_{1}=m\geq 0$, ${S}_{A}=({S}_{A_{1}}^{\prime},\ldots,{S}_{A_{m}}^{\prime})^{\prime}$, ${S}_{A_{k}}=\lambda\sqrt{d_{A_{k}}}U_{A_{k}}$ and $\|U_{A_{k}}\|_{2}=1$.

For a sufficiently large constant $C_{2}>0$, define

As in the proof of Theorem 2.1 of Wei and Huang (2008),

for a constant $M_{1}>1$. By the triangle and Cauchy–Schwarz inequalities,

In the proof of Theorem 2.1 of Wei and Huang (2008), it is shown that

and $m_{n}=O(n^{1/(2d+1)})$, we have for all $t\geq 0$ and $n$ sufficiently large,

Before proving part (ii), we first prove part (iii) of Theorem 1. By the definition of $\widetilde{\bolds\beta}_{n}\equiv(\widetilde{\bolds\beta}_{n1}^{\prime},\ldots,\widetilde{\bolds\beta}_{np}^{\prime})^{\prime}$,

Let $A_{2}=\{j\dvtx\|\bolds\beta_{nj}\|_{2}\neq 0\mbox{ or }\|\widetilde{\bolds\beta}_{nj}\|_{2}\neq 0\}$ and $d_{n2}=|A_{2}|$. By part (i), $d_{n2}=O_{p}(q)$. By (19) and the definition of $A_{2}$,

Let $\bolds\eta_{n}={\mathbf{Y}}-{\mathbf{Z}}\bolds\beta_{n}$. Write

Let $\bolds\nu_{n}={\mathbf{Z}}_{A_{2}}(\widetilde{\bolds\beta}_{nA_{2}}-\bolds\beta_{nA_{2}})$. Combining (Appendix: Proofs), (Appendix: Proofs) and (Appendix: Proofs) to get

Let $\bolds\eta_{n}^{*}$ be the projection of $\bolds\eta_{n}$ to the span of ${\mathbf{Z}}_{A_{2}}$, that is, $\bolds\eta_{n}^{*}={\mathbf{Z}}_{A_{2}}({\mathbf{Z}}_{A_{2}}^{\prime}\times\break{\mathbf{Z}}_{A_{2}})^{-1}{\mathbf{Z}}_{A_{2}}^{\prime}\bolds\eta_{n}$. By the Cauchy–Schwarz inequality,

Let $c_{n*}$ be the smallest eigenvalue of ${\mathbf{Z}}_{A_{2}}^{\prime}{\mathbf{Z}}_{A_{2}}/n$. By Lemma 3 and part (i), $c_{n*}\asymp_{p}m_{n}^{-1}$. Since $\|\bolds\nu_{n}\|_{2}^{2}\geq nc_{n*}\|\widetilde{\bolds\beta}_{nA_{2}}-\bolds\beta_{nA_{2}}\|_{2}^{2}$ and $2ab\leq a^{2}+b^{2}$,

Let $f_{0}({\mathbf{X}}_{i})=\sum_{j=1}^{p}f_{0j}(X_{ij})$ and $f_{0A}({\mathbf{X}}_{i})=\sum_{j\in A}f_{0j}(X_{ij})$. Write

Since $|\mu-\overline{Y}|^{2}=O_{p}(n^{-1})$ and $\|f_{0j}-f_{nj}\|_{\infty}=O(m_{n}^{-d})$, we have

where $\bolds\varepsilon_{n}^{*}$ is the projection of $\bolds\varepsilon_{n}=(\varepsilon_{1},\ldots,\varepsilon_{n})^{\prime}$ to the span of ${\mathbf{Z}}_{A_{2}}$. We have

where $\mathcal{Z}_{jk}=(\psi_{k}(X_{1j}),\ldots,\psi_{k}(X_{nj}))^{\prime}$. By Lemma 2,

Since $d_{n2}=O_{p}(q)$, $c_{n*}\asymp_{p}m_{n}^{-1}$ and $c_{n}^{*}\asymp_{p}m_{n}^{-1}$, we have

We now prove part (ii). Since $\|f_{j}\|_{2}\geq c_{f}>0,1\leq j\leq q$, $\|f_{j}-f_{nj}\|_{2}=O(m_{n}^{-d})$ and $\|f_{nj}\|_{2}\geq\|f_{j}\|_{2}-\|f_{j}-f_{nj}\|_{2}$, we have $\|f_{nj}\|_{2}\geq 0.5c_{f}$ for $n$ sufficiently large. By a result of de Boor (2001), see also (12) of Stone (1986), there are positive constants $c_{6}$ and $c_{7}$ such that

It follows that $\|\bolds\beta_{nj}\|_{2}^{2}\geq c_{7}^{-1}m_{n}\|f_{nj}\|_{2}^{2}\geq 0.25c_{7}^{-1}c_{f}^{2}m_{n}.$ Therefore, if $\|\bolds\beta_{nj}\|_{2}\neq 0$ but $\|\widetilde{\bolds\beta}_{nj}\|_{2}=0$, then

However, since $(m_{n}\log(pm_{n}))/n\rightarrow 0$ and $(\lambda_{n1}^{2}m_{n})/n^{2}\rightarrow$, (28) contradicts part (iii). {pf*}Proof of Theorem 2 By the definition of $\widetilde{f}_{j},1\leq j\leq p$, parts (i) and (ii) follow from parts (i) and (ii) of Theorem 1 directly.

Now consider part (iii). By the properties of spline [de Boor (2001)],

Part (iii) follows from (Appendix: Proofs) and (30).

In the proofs below, for any matrix ${\mathbf{H}}$, denote its $2$-norm by $\|{\mathbf{H}}\|$, which is equal to its largest eigenvalue. This norm satisfies the inequality $\|{\mathbf{H}}{\mathbf{x}}\|\leq\|{\mathbf{H}}\|\|{\mathbf{x}}\|$ for a column vector ${\mathbf{x}}$ whose dimension is the same as the number of the columns of ${\mathbf{H}}$.

Denote $\bolds\beta_{nA_{1}}=(\bolds\beta_{nj}^{\prime},j\in A_{1})^{\prime}$, $\widehat{\bolds\beta}_{nA_{1}}=(\widehat{\bolds\beta}_{nj}^{\prime},j\in A_{1})^{\prime}$ and ${\mathbf{Z}}_{A_{1}}=({\mathbf{Z}}_{j},j\in A_{1})$. Define ${\mathbf{C}}_{A_{1}}=n^{-1}{\mathbf{Z}}_{A_{1}}^{\prime}{\mathbf{Z}}_{A_{1}}$. Let $\rho_{n1}$ and $\rho_{n2}$ be the smallest and largest eigenvalues of ${\mathbf{C}}_{A_{1}}$, respectively. {pf*}Proof of Theorem 3 By the KKT, a necessary and sufficient condition for $\widehat{\bolds\beta}_{n}$ is

Let $\bolds\nu_{n}=(w_{nj}\widehat{\bolds\beta}_{j}/(2\|\widehat{\bolds\beta}_{nj}\|),j\in A_{1})^{\prime}$. Define

If $\widehat{\bolds\beta}_{nA_{1}}=_{0}\bolds\beta_{nA_{1}}$, then the equation in (31) holds for $\widehat{\bolds\beta}_{n}\equiv(\widehat{\bolds\beta}_{nA_{1}}^{\prime},{\mathbf{0}}^{\prime})^{\prime}$. Thus, since ${\mathbf{Z}}\widehat{\bolds\beta}_{n}={\mathbf{Z}}_{A_{1}}\widehat{\bolds\beta}_{nA_{1}}$ for this $\widehat{\bolds\beta}_{n}$ and $\{{\mathbf{Z}}_{j},j\in A_{1}\}$ are linearly independent,

Let $f_{0j}({\mathbf{X}}_{j})=(f_{0j}(X_{1j}),\ldots,f_{0j}(X_{nj}))^{\prime}$ and $\bolds\delta_{n}=\sum_{j\in A_{1}}f_{0j}({\mathbf{X}}_{j})-\break{\mathbf{Z}}_{A_{1}}\bolds\beta_{nA_{1}}$. By Lemma 1, we have

Let ${\mathbf{H}}_{n}={\mathbf{I}}_{n}-{\mathbf{Z}}_{A_{1}}({\mathbf{Z}}_{A_{1}}^{\prime}{\mathbf{Z}}_{A_{1}})^{-1}{\mathbf{Z}}_{A_{1}}^{\prime}$. By (32),

Based on these two equations, Lemma 5 below shows that

These two equations lead to part (i) of the theorem.

We now prove part (ii) of Theorem 3. As in (26), for $\bolds\eta_{n}={\mathbf{Y}}-{\mathbf{Z}}\bolds\beta_{n}$ and

where $\bolds\varepsilon_{n1}^{*}$ is the projection of $\bolds\varepsilon_{n}=(\varepsilon_{1},\ldots,\varepsilon_{n})^{\prime}$ to the span of ${\mathbf{Z}}_{A_{1}}$. We have

Now similarly to the proof of (25), we can show that

Since $\rho_{n1}\asymp_{p}m_{n}^{-1}$, the result follows.

The following lemmas are needed in the proof of Theorem 3.

For $\bolds\nu_{n}=(w_{nj}\widetilde{\bolds\beta}_{j}/(2\|\widetilde{\bolds\beta}_{nj}\|),j\in A_{1})^{\prime}$, under condition (B1),

and $\sum_{j\in A_{1}}\|\bolds\beta_{nj}\|^{-2}\leq qb_{n1}^{-2}.$ The claim follows.

Let $\rho_{n3}$ be the maximum of the largest eigenvalues of $n^{-1}{\mathbf{Z}}_{j}^{\prime}{\mathbf{Z}}_{j},j\in A_{0}$, that is, $\rho_{n3}=\max_{j\in A_{0}}\|n^{-1}{\mathbf{Z}}_{j}^{\prime}{\mathbf{Z}}_{j}\|_{2}$. By Lemma 3,

Under conditions (B1), (B2), (A3) and (A4),

Let ${\mathbf{T}}_{nj}$ be an $m_{n}\times qm_{n}$ matrix with the form

where ${\mathbf{O}}_{m_{n}}$ is an $m_{n}\times m_{n}$ matrix of zeros and ${\mathbf{I}}_{m_{n}}$ is an $m_{n}\times m_{n}$ identity matrix, and ${\mathbf{I}}_{m_{n}}$ is at the $j$th block. By (34), $\widehat{\bolds\beta}_{nj}-\bolds\beta_{nj}=n^{-1}{\mathbf{T}}_{nj}{\mathbf{C}}_{A_{1}}^{-1}({\mathbf{Z}}_{A_{1}}^{\prime}\bolds\varepsilon_{n}+{\mathbf{Z}}_{A_{1}}^{\prime}\bolds\delta_{n}-\lambda_{n2}\bolds\nu_{n}).$ By the triangle inequality,

Let $C$ be a generic constant independent of $n$. The first term on the right-hand side

Thus, (40) follows from (39), (Appendix: Proofs)–(44) and condition (B2a).

Under conditions (B1), (B2), (A3) and (A4),

Recall $s_{n}=p-q$ is the number of zero components in the model. By Lemma 2,

Since $w_{nj}=\|\widehat{\bolds\beta}_{nj}\|^{-1}=O_{p}(r_{n})$ for $j\notin A_{1}$ and by (47), for the first term on the right-hand side of (46), we have

By (33), the second term on the right-hand side of (46)

By Lemma 4, the third term on the right-hand side of (46)

Therefore, (45) follows from (39), (Appendix: Proofs), (Appendix: Proofs), (Appendix: Proofs) and condition (B2b). {pf*}Proof of Theorem 4 The proof is similar to that of Theorem 2 and is omitted.

Acknowledgments

The authors wish to thank the Editor, Associate Editor and two anonymous referees for their helpful comments.

References