Prediction Weighted Maximum Frequency Selection

TL;DR

This paper introduces a novel variable selection method combining adaptive LASSO and Elastic-Net estimators with a prediction-weighted maximum frequency approach, achieving consistent model selection and improved finite sample performance.

Contribution

It develops a new strategy for variable selection that enhances finite sample performance and consistency using prediction-weighted maximum frequency models with adaptive estimators.

Findings

Achieves consistent model selection with divergence of p_n.

Demonstrates excellent finite sample performance in simulations.

Extends methodology to generalized linear models.

Abstract

Shrinkage estimators that possess the ability to produce sparse solutions have become increasingly important to the analysis of today's complex datasets. Examples include the LASSO, the Elastic-Net and their adaptive counterparts. Estimation of penalty parameters still presents difficulties however. While variable selection consistent procedures have been developed, their finite sample performance can often be less than satisfactory. We develop a new strategy for variable selection using the adaptive LASSO and adaptive Elastic-Net estimators with $p_n$ diverging. The basic idea first involves using the trace paths of their LARS solutions to bootstrap estimates of maximum frequency (MF) models conditioned on dimension. Conditioning on dimension effectively mitigates overfitting, however to deal with underfitting, these MFs are then prediction-weighted, and it is shown that not only can consistent model selection be achieved, but that attractive convergence rates can as well, leading to excellent finite sample performance. Detailed numerical studies are carried out on both simulated and real datasets. Extensions to the class of generalized linear models are also detailed.

Tables

Table 1. Table 1: The leukaemia classification using adaptive LASSO

Criteria	Ten-fold CV error	Test error	Number of genes
WMF	0/38	5/34	5
CV	0/38	4/34	13
Cp	0/38	4/34	18
BIC	0/38	4/34	18
EBIC	1/38	6/34	5
GIC	1/18	6/34	5

Table 2. Table 2: The leukaemia classification using adaptive Elastic-Net

Criteria	Ten-fold CV error	Test error	Number of genes
WMF	0/38	4/34	10
CV	1/38	6/34	42
Cp	1/38	6/34	36
BIC	1/38	7/34	34
EBIC	1/38	7/34	34
GIC	1/38	7/34	21

Equations

$$\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{\epsilon },$$

$$\hat{\mathbf{\beta }}=\arg \underset{\mathbf{\beta }}{\min }\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }^2+\sum _{j=1}^{p_n}\rho (|{\beta }_j|,\mathbf{\lambda }),$$

$${\hat{\mathbf{\beta }}}_a=\arg \underset{\mathbf{\beta }}{\min }||\mathbf{y}-\mathbf{X}\mathbf{\beta }|{|}^2+2{\lambda }_n\sum _{j=1}^{p_n}{\omega }_j|{\beta }_j|,$$

$${\hat{\mathbf{\beta }}}_e=(1+\frac{{\lambda }_{n2}}{n})\left\{\arg \underset{\mathbf{\beta }}{\min }||\mathbf{y}-\mathbf{X}\mathbf{\beta }|{|}^2+{\lambda }_{n2}\sum _{j=1}^{p_n}|{\beta }_j{|}^2+{\lambda }_{n1}\sum _{j=1}^{p_n}|{\beta }_j|\right\}.$$

$${\hat{\mathbf{\beta }}}_{ae}=(1+\frac{{\lambda }_{n2}}{n})\left\{\arg \underset{\mathbf{\beta }}{\min }||\mathbf{y}-\mathbf{X}\mathbf{\beta }|{|}^2+{\lambda }_{n2}\sum _{j=1}^{p_n}|{\beta }_j{|}^2+{\lambda }_{n1}^{+}\sum _{j=1}^{p_n}{\omega }_j|{\beta }_j|\right\},$$

$${\hat{\epsilon }}_{0i}=y_i-{\mathbf{x}}_i^T\hat{\mathbf{\beta }},$$

$$r^{\ast }=\max \{j:j=\arg \underset{1\le i\le p-1}{\max }{\text{MF}}_i\}.$$

$${\hat{\beta }}_j=\{|{\mathbf{X}}_j^T\mathbf{y}|-\frac{{\lambda }_n}{|{\tilde{\beta }}_j{|}^{\gamma }}{\}}_{+}sgn({\mathbf{X}}_j^T\mathbf{y}),j=1,\dots ,p_n,$$

$${\mathbf{X}}_j^T\mathbf{y}={\beta }_{0j}+{\mathbf{X}}_j^T\mathbf{\epsilon },$$

$$P\left(|{\mathbf{X}}_i^T\mathbf{y}|>|{\mathbf{X}}_j^T\mathbf{y}|\right)>0.5\text{if}|{\beta }_{0i}|>|{\beta }_{0j}|,$$

$$P\left(|{\mathbf{X}}_i^T\mathbf{y}|>|{\mathbf{X}}_j^T\mathbf{y}|\right)=0.5\text{if}|{\beta }_{0i}|=|{\beta }_{0j}|,$$

$$(1).P({\hat{\mathcal{A}}}_r={\mathcal{S}}_r\mid r)>P({\hat{\mathcal{A}}}_r={\mathcal{W}}_r\mid r),0<r\le p_0,$$

$$(2).P({\hat{\mathcal{A}}}_r={\mathcal{W}}_r^1\mid r)=P({\hat{\mathcal{A}}}_r={\mathcal{W}}_r^2\mid r),p_0<r<p_n,$$

$$d\le {\zeta }_{min}(\frac{1}{n}{\mathbf{X}}^T\mathbf{X})\le {\zeta }_{max}(\frac{1}{n}{\mathbf{X}}^T\mathbf{X})\le D,$$

$$\underset{n\to \infty }{\lim }\frac{p_n}{n}\frac{1}{{\min }_{j\in \mathcal{A}}|{\beta }_{0j}{|}^2}\to 0.$$

$$\underset{n\to \infty }{\lim }{\lambda }_n/\sqrt{n}\to 0,\underset{n\to \infty }{\lim }\frac{{\lambda }_n}{\sqrt{n}}{n}^{\frac{(1-\varrho )(1+\gamma )-1}{2}}\to \infty ,$$

$$\underset{n\to \infty }{\lim }{\left(\frac{{\lambda }_n}{\sqrt{n}}\right)}^{\frac{1}{\gamma }}\frac{1}{{\min }_{j\in \mathcal{A}}|{\beta }_{0j}|}\to 0.$$

$$\underset{n\to \infty }{\lim }{\lambda }_{n1}/\sqrt{n}\to 0,\underset{n\to \infty }{\lim }{\lambda }_{n2}/\sqrt{n}\to 0,$$

$$\underset{n\to \infty }{\lim }{\lambda }_{n1}^{+}/\sqrt{n}\to 0,\underset{n\to \infty }{\lim }\frac{{\lambda }_{n1}^{+}}{\sqrt{n}}{n}^{\frac{(1-\varrho )(1+\gamma )-1}{2}}\to \infty ,$$

$$\underset{n\to \infty }{\lim }{\left(\frac{{\lambda }_{n_1}^{+}}{\sqrt{n}}\right)}^{\frac{1}{\gamma }}\frac{1}{{\min }_{j\in \mathcal{A}}|{\beta }_{0j}|}\to 0.$$

$$\underset{n\to \infty }{\lim }{P}^{\ast }({\mathcal{A}}_n^{\ast }=\mathcal{A}\mid {\lambda }_n)=1.$$

$$\underset{n\to \infty }{\lim }{P}^{\ast }({\mathcal{A}}_n^{\ast }={\mathcal{M}}_r\mid {\lambda }_n^{\prime })<1,$$

$${\lambda }_n(0)>{\lambda }_n(1)>\cdots >{\lambda }_n(10)=0.$$

$${\lambda }_n(m+1)=\arg \underset{{\lambda }_n}{\min }\Vert \mathbf{y}-\mathbf{X}{\hat{\mathbf{\beta }}}_a({\lambda }_n){\Vert }^2+a_n\widehat{df}({\lambda }_n),{\lambda }_n(m+1)\le {\lambda }_n<{\lambda }_n(m),$$

$$\Vert \mathbf{y}-\mathbf{X}{\hat{\mathbf{\beta }}}_a({\lambda }_n(m_k)){\Vert }^2<\Vert \mathbf{y}-\mathbf{X}{\hat{\mathbf{\beta }}}_a({\lambda }_n(m_k^{\prime })){\Vert }^2.$$

$$k⟺{\lambda }_n(m_k).$$

$$\underset{n\to \infty }{\lim }{P}^{\ast }({\mathcal{A}}_n^{\ast }=\mathcal{A}\mid p_0)=1,$$

$$\underset{n\to \infty }{\lim }{P}^{\ast }({\mathcal{A}}_n^{\ast }={\mathcal{M}}_r\mid r)<1,p_0<r<p_n,$$

$$\underset{n\to \infty }{\lim }{P}^{\ast }({\mathcal{T}}_n^{\ast }=\mathcal{A}\mid p_0)=1,$$

$$\underset{n\to \infty }{\lim }{P}^{\ast }({\mathcal{T}}_n^{\ast }={\mathcal{M}}_r\mid r)<1,p_0<r<p_n,$$

$$\underset{n\to \infty }{\lim }{P}^{\ast }({\mathcal{T}}_n^{\ast }=\mathcal{A}\mid {\lambda }_{n1}^{+})=1.$$

$$\underset{n\to \infty }{\lim }{P}^{\ast }({\mathcal{T}}_n^{\ast }={\mathcal{M}}_r\mid {\lambda }_{n1}^{\prime })<1,$$

Full text

**Prediction Weighted Maximum Frequency Selection

**Hongmei Liu and J. Sunil Rao *Division of Biostatistics, University of Miami, Miami, FL 33136

h.liu7@med.miami.edu

JRao@biostat.med.miami.edu*

Shrinkage estimators that possess the ability to produce sparse solutions have become increasingly important to the analysis of today’s complex datasets. Examples include the LASSO, the Elastic-Net and their adaptive counterparts. Estimation of penalty parameters still presents difficulties however. While variable selection consistent procedures have been developed, their finite sample performance can often be less than satisfactory. We develop a new strategy for variable selection using the adaptive LASSO and adaptive Elastic-Net estimators with $p_{n}$ diverging. The basic idea first involves using the trace paths of their LARS solutions to bootstrap estimates of maximum frequency (MF) models conditioned on dimension. Conditioning on dimension effectively mitigates overfitting, however to deal with underfitting, these MFs are then prediction-weighted, and it is shown that not only can consistent model selection be achieved, but that attractive convergence rates can as well, leading to excellent finite sample performance. Detailed numerical studies are carried out on both simulated and real datasets. Extensions to the class of generalized linear models are also detailed. Key Words: Adaptive LASSO, adaptive Elastic-Net, model selection, bootstrap.

1 Introduction

Consider the standard linear regression model

[TABLE]

where $\mathbf{y}=(y_{1},\dots,y_{n})^{T}$ is a vector of responses, $\mathbf{X}=(\mathbf{X}_{1},\dots,\mathbf{X}_{p_{n}})$ is an $n\times p_{n}$ design matrix of predictors, ${\boldsymbol{\beta}}=(\beta_{1},\dots,\beta_{p_{n}})^{T}$ is a vector of unknown regression parameters, $\boldsymbol{\varepsilon}=(\varepsilon_{1},\dots,\varepsilon_{n})^{T}$ is a vector of independent and identically distributed (i.i.d.) random errors. We allow $p_{n}$ to increase with $n$.

Because some elements of ${\boldsymbol{\beta}}$ might be [math], a family of penalized least squares estimators were developed for variable selection and estimation,

[TABLE]

where $\|\cdot\|$ is the $L_{2}$-norm, $\boldsymbol{\lambda}\geq\mathbf{0}$ are regularization parameters, and $\rho(|\beta_{j}|,\boldsymbol{\lambda})$ is positive valued for $\beta_{j}\neq 0$. [1] pointed out that through variable selection one can focus on a small number of important predictors for enhanced scientific discovery and potentially improve prediction performance by removing noise variables.

Penalized $L_{q}$-regression is a special case of (1.2) with $\rho(|\beta_{j}|,\boldsymbol{\lambda})=\lambda|\beta_{j}|^{q},q\geq 0$, which includes the best subset selection for $q=0$; the LASSO [2] for $q=1$ and the ridge regression [3] for $q=2$. Best subset selection is known to be computationally infeasible for high dimensional data and inherently discrete in variable selection [4]. Ridge regression does not possess a variable selection property. The LASSO however, can do simultaneous estimation and variable selection because its $L_{1}$ penalty is singular at the origin and can shrink some coefficients to exact 0 with a sufficiently large $\lambda$ [5]. Other penalized least squares estimators that can do simultaneous estimation and variable selection include the SCAD [5] and adaptive LASSO [6] both enjoying the oracle properties [5]; the Elastic-Net [7] capable of detecting grouped effects; the adaptive Elastic-Net [8] combining advantages of the adaptive LASSO and Elastic-Net; and etc.

Selection of $\boldsymbol{\lambda}$ is essential in above penalized least squares estimation procedures. Although methods such as the SCAD, adaptive LASSO and adaptive Elastic-Net enjoy the oracle properties asymptotically, their optimal properties rely on particular specifications of the $\boldsymbol{\lambda}$, whose magnitude controls the complexity of a selected model and trade-off between bias and variance in estimators [9]. The multi-fold cross-validation (CV) and generalized cross-validation (GCV) are frequently applied for the tuning parameters selection [2, 5, 6, 7]. But they overfit the model asymptotically [10]. For consistent variable selection, [11] suggested to use the BIC in adaptive LASSO and a modified BIC when $p_{n}$ is diverging [12]; [13] introduced an extended BIC (EBIC) for linear models and then generalized it to generalized linear models [14]; [15] put forward a generalized information criterion (GIC) with $p_{n}$ diverging; [16] provided a consistent cross-validation procedure (CCV) for the LASSO; [17] proposed the stability selection (SS) for their randomized LASSO. Although variable selection consistency was established for these procedures, their finite sample performance can often be less than optimal (Section 6 ahead demonstrates this in simulation studies).

We propose a new method for tuning parameters selection, focusing in particular on the adaptive LASSO and adaptive Elastic-Net estimators. A simple example helps to illustrate the basic idea. Consider the adaptive LASSO in following example.

Example 1. Data are drawn from model (1.1) with $\boldsymbol{\beta}=(3,1.5,0,0,2,0,\dots)_{10}^{T}$, row vectors of the design matrix $\mathbf{x}_{i}\overset{iid}{\sim}N_{10}(\mathbf{0},\boldsymbol{\Sigma})$ with $\boldsymbol{\Sigma}(i,j)=0.3^{|i-j|}$ and $\varepsilon_{i}\overset{iid}{\sim}N(0,3^{2})$ for $i=1,\dots,100$. So the true model size here is 3.

Figure 1 (top) shows the full adaptive LASSO solution path from the LARS algorithm [18]. In the figure, each step indicates a dimension change in the estimator. These steps are called transition points in [19]. They showed that if using information criteria such as the AIC or BIC to identify the optimal $\lambda$ in adaptive LASSO, it lies in one of the transition points. This result helps to justify uses of the LARS algorithm and our subsequent focus on the transition points. Then the question remains about how to choose from these transition points.

Figure 1 (middle and bottom) gives a brief look at our proposed method. The middle panel shows the estimated maximum frequency (MF) of a candidate model given the dimension. The MF estimation is done by a bootstrapping algorithm using the transition points. The strategy of conditioning on dimension has two important consequences: i) for overfit dimensions, the MFs are dramatically smaller than the true dimension MF (other than the full model of course), and ii) underfit dimensions can also produce large MF values. Point i) is important because for variable selection, overfitting is usually much more difficult to deal with. So one must now deal with the underfitting issue. We do this by introducing a prediction-based weight to the MFs (labeled as WMF). The results are shown in bottom panel of Figure 1. As is evident, now the true dimension, which maps to the true model, stands out beautifully from all others.

The rest of the paper is organized as follows. In Section 2, we briefly review the adaptive LASSO and adaptive Elastic-Net estimators and introduce the bootstrap algorithm for each. In Section 3, the MF procedure itself is described and its underlying properties are carefully examined using a simple orthogonal design. In Section 4, asymptotic properties of the MF procedure are established in general settings. The WMF procedure and its variable selection consistency are presented in Section 5. Comprehensive simulation studies are shown in Section 6. Section 7 describes extensions of the MWF procedure to generalized linear models (GLMs). Applications of the WMF procedure to ultra-high dimensional data are discussed in Section 8.

2 Bootstrapping the adaptive LASSO and adaptive Elastic-Net estimators

Denote $\boldsymbol{\beta}_{0}$ the true value of $\boldsymbol{\beta}$ with model size $p_{0}$, and $\boldsymbol{\tilde{\beta}}=(\tilde{\beta}_{1},\dots,\tilde{\beta}_{p_{n}})^{T}$ a consistent estimate of $\boldsymbol{\beta}_{0}$. The adaptive LASSO [6] estimator is

[TABLE]

where $\omega_{j}=|\tilde{\beta}_{j}|^{-\gamma}$, $\gamma\geq 0$. It was suggested to use the ordinary least-squares (OLS) estimator or the best ridge estimator (if collinearity exists) for $\boldsymbol{\tilde{\beta}}$. Under certain regularity conditions, $\boldsymbol{\hat{\beta}}_{a}$ was shown to enjoy the oracle properties.

The Elastic-Net estimator [7] is

[TABLE]

It overcomes several limitations pertaining to the LASSO: 1) the added $L_{2}$ penalty is strictly convex to allow grouping effects; 2) In a $p_{n}>n$ case, it can potentially estimate all $p_{n}$ predictors, while the LASSO can only find at most $n$ predictors.

[8] proposed the adaptive Elastic-Net to combine strengths of the Elastic-Net and adaptive LASSO. The adaptive Elastic-Net estimator is

[TABLE]

where $\omega_{j}=|\hat{\beta}_{ej}|^{-\gamma}$, $\gamma\geq 0$ and $\boldsymbol{\hat{\beta}}_{e}=(\hat{\beta}_{e1},\dots,\hat{\beta}_{ep_{n}})^{T}$ is the Elastic-Net estimator in (2.2). Note that $\lambda_{n2}$ takes the same value for the $L_{2}$ penalty function in (2.2) and (2.3), because the $L_{2}$ penalty contributes to the same kind of grouping effects. On the other hand, $\lambda_{n1}$ and $\lambda_{n1}^{+}$ are allowed to be different as they control the sparsity in estimators. Under some regularity conditions, $\boldsymbol{\hat{\beta}}_{ae}$ was shown to enjoy the oracle properties.

We now detail bootstrapping for these two estimators. There are typically two ways of generating bootstrap observations for model (1.1) [20].

1. Bootstrapping pairs [21]. Let $\hat{F}(\mathbf{X},\mathbf{y})$ be the empirical distribution putting mass $n^{-1}$ on each data pair $(\mathbf{x}_{i},y_{i}),i=1,\dots,n$. Generate i.i.d. paired bootstrap data $\{(\mathbf{x}_{i}^{*},y_{i}^{*}),i=1,\dots,n\}$ from $\hat{F}(\mathbf{X},\mathbf{y})$. The bootstrap analog of $\boldsymbol{\hat{\beta}}_{a}$, denoted as $\boldsymbol{\hat{\beta}}_{a}^{*}$, is to replace $(\mathbf{X},\mathbf{y})$ with $(\mathbf{X}^{*},\mathbf{y}^{*})$ in (2.1) where $\mathbf{X}^{*}=(\mathbf{x}_{1}^{*},\dots,\mathbf{x}_{n}^{*})^{T}$ and $\mathbf{y}^{*}=(y_{1}^{*},\dots,y_{n}^{*})^{T}$. So is the bootstrap analog of $\boldsymbol{\hat{\beta}}_{ae}$, denoted as $\boldsymbol{\hat{\beta}}_{ae}^{*}$. Under the weak condition that $\mathbf{X}^{T}\mathbf{X}\to\infty$, $\mathbf{X}^{T}\mathbf{X}(\mathbf{X}^{*T}\mathbf{X}^{*})^{-1}\to 1$ almost surely [20].

2. Bootstrapping residuals [22]. Calculate the $i$th residual

[TABLE]

where $\boldsymbol{\hat{\beta}}$ is a ridge estimate of $\boldsymbol{\beta}_{0}$. Generate i.i.d. bootstrap residuals $\{\varepsilon_{i}^{*},i=1,\dots,n\}$ from the empirical distribution that puts mass $n^{-1}$ on each centered residual, $\hat{\varepsilon}_{i}=\hat{\varepsilon}_{0i}-\bar{\varepsilon}_{0}$, where $\bar{\varepsilon}_{0}$ is the average of $\hat{\varepsilon}_{0i},i=1,\dots,n$. Then the i.i.d. residual bootstrap data is $\{(\mathbf{x}_{i},y_{i}^{*}),i=1,\dots,n\}$ where $y_{i}^{*}=\mathbf{x}_{i}^{T}\boldsymbol{\hat{\beta}}+\varepsilon_{i}^{*}$. The bootstrap analog of $\boldsymbol{\hat{\beta}}_{a}$, denoted as $\boldsymbol{\tilde{\beta}}_{a}^{*}$, is to substitute $\mathbf{y}$ with $\mathbf{y}^{*}$ in (2.1). So is the bootstrap analog of $\boldsymbol{\hat{\beta}}_{ae}$, denoted as $\boldsymbol{\tilde{\beta}}_{ae}^{*}$.

In next section, we introduce the MF procedure which takes use of above bootstrap estimators.

3 The MF procedure

Denote a $j$-dimensional candidate model from the $i$th bootstrap data as $M_{j}^{i}$.

Remark 1**.**

In the 4th step, the full model is excluded because it will destroy the maximum frequency rule by having the highest frequency, $B$, all the time. If there is a tie at the maximum of $\text{MF}_{i},1\leq i\leq p-1$, we select the one at the highest dimension. This strategy guarantees asymptotic variable selection consistency of the MF procedure, which will be discussed in Section 4.

The MF procedure for adaptive Elastic-Net is in parallel. But in the 2nd step, the LARS-EN algorithm [7] is used instead to fit each bootstrap data.

We discussed in introduction to this paper consequences of the MF procedure by conditioning on dimension. Here we use a simple orthogonal design with i.i.d. normal random errors to study underlying properties driving that performance. In this case, we have $\mathbf{X}^{T}\mathbf{X}=\mathbf{I}$ and the adaptive Elastic-Net reduces automatically to the adaptive LASSO [8]. Denote $\mathbf{X}_{j}$ the $j$th column of $\mathbf{X}$. Then the adaptive LASSO estimator is

[TABLE]

where $z_{+}$ equals to $z$ if $z>0$ otherwise 0. We can expand the $\mathbf{X}_{j}^{T}\mathbf{y}$ by

[TABLE]

where $\mathbf{X}_{j}^{T}{\boldsymbol{\varepsilon}}\sim N(0,\sigma^{2})$. The following Lemma gives an order relationship for $\mathbf{X}_{j}^{T}\mathbf{y}$’s.

Lemma 1**.**

Suppose ${\mathbf{X}}^{T}{\mathbf{X}}={\mathbf{I}}$, then we have

[TABLE]

for $i,j\in\{1,\dots,p_{n}\}$.

In combine with the fact that $\frac{\lambda_{n}}{|\tilde{\beta}_{i}|^{\gamma}}>\frac{\lambda_{n}}{|\tilde{\beta}_{j}|^{\gamma}}$ asymptotically for $\beta_{0i}<\beta_{0j}$, it is easy to deduce from (3.1) that given a $\lambda_{n}$ adaptive LASSO tends to select those variables, corresponding to the first $k_{\lambda_{n}}$ largest $|\beta_{j}|$’s, with the highest probability.

Without loss of generality, suppose $|\boldsymbol{\beta}_{0}|$ is decreasingly ordered. Denote $\mathcal{S}_{r}$ a $r$-dimensional model containing the first $r$ elements of $|\boldsymbol{\beta}_{0}|$, and denote $\mathcal{W}_{r}$ any other $r$-dimensional models. Let $\hat{\mathcal{A}}_{r}$ be an adaptive LASSO model estimate given the model size is $r$, $P(\hat{\mathcal{A}}_{r}=\mathcal{S}_{r}\mid r)$ indicates the conditional probability of $\hat{\mathcal{A}}_{r}=\mathcal{S}_{r}$ given the model size. Then preceding deductions from (3.1) can be formularized as

[TABLE]

where $\mathcal{W}_{r}^{1}$ and $\mathcal{W}_{r}^{2}$ are two $r$-dimensional models s.t. $\mathcal{S}_{p_{0}}\subset\mathcal{W}_{r}^{1},\mathcal{W}_{r}^{2}$.

Above properties of the adaptive LASSO coincides to some extent with the results of Theorem 2 in [23]. By (3.3), zero predictors will be equally likely selected at an overfit dimension. As a result $P(\hat{\mathcal{A}}_{r}=M_{r}\mid r)$ (see Algorithm 1 for definition of $M_{r}$), $p_{0}<r<p_{n}$, drops down dramatically, which is why we see a huge gap between the true dimension and overfit dimensions in Figure 1 (middle). On the other hand, $P(\hat{\mathcal{A}}_{r}=\mathcal{S}_{r}\mid r)$ at some underfit dimensions can be as competitive as $P(\hat{\mathcal{A}}_{r}=\mathcal{S}_{p_{0}}\mid p_{0})$. We propose a WMF procedure to tackle this underfitting issue in Section 5.

In next section, we show asymptotic variable selection properties for $\boldsymbol{\hat{\beta}}_{a}^{*}$ and $\boldsymbol{\tilde{\beta}}_{ae}^{*}$ in general settings, from which variable selection consistency of the MF procedure can be deduced.

4 Asymptotic properties of the MF procedure

Let $\mathcal{A}=\{j:\beta_{0j}\neq 0\}$ be the true model. We assume following regularity conditions for subsequent theoretical studies:

(A1) Denote $\zeta_{min}(\mathbf{C})$ and $\zeta_{max}(\mathbf{C})$ the minimum and maximum eigenvalues of a positive definite matrix $\mathbf{C}$. We assume

[TABLE]

where $d$ and $D$ are two positive constants.

(A2) $p_{n}=n^{\varrho},\ 0\leq\varrho<1$ and $\gamma>{\varrho\over 1-\varrho}$. The last inequation is to ensure $(1-\varrho)(1+\gamma)>1$ in (A3)–(A4). Moreover,

[TABLE]

(A3) In adaptive LASSO,

[TABLE]

and

[TABLE]

(A4) In adaptive Elastic-Net,

[TABLE]

and

[TABLE]

[TABLE]

(A5) The errors $\{\varepsilon_{i},i=1,\dots,n\}$ are i.i.d. with mean 0 and variance $\sigma^{2}<\infty$.

Denote $\mathcal{A}_{n}^{*}=\{j:\hat{\beta}_{aj}^{*}\neq 0\}$ an adaptive LASSO estimate of $\mathcal{A}$ using paired bootstrap data. Let $P^{*}=P(\cdot\mid\mathcal{E})$ and $E^{*}=E(\cdot\mid\mathcal{E})$ where $\mathcal{E}=\sigma\left((\mathbf{x}_{i},y_{i}),i=1,\dots,n\right)$. Then $P^{*}(\mathcal{A}_{n}^{*}=\mathcal{A}\mid\lambda_{n})$ indicates the conditional probability of $\mathcal{A}_{n}^{*}=\mathcal{A}$ given $\mathcal{E}$ and $\lambda_{n}$.

Theorem 1**.**

Suppose conditions (A1)–(A3) and (A5) hold, then

[TABLE]

Moreover, let $\lambda_{n}^{\prime}$ be another tuning parameter such that the adaptive LASSO estimator under $\lambda_{n}^{\prime}$ is of dimension $r$, $p_{0}<r<p_{n}$, then

[TABLE]

where $\mathcal{M}_{r}$ is any $r$-dimensional model.

Proofs of Theorem 1 are included in Appendix A.

In adaptive LASSO, given a $\lambda_{n}$ is equivalent to given a dimension, but the converse is not true. One dimension can be mapped to numerous models, as a result to numerous tuning parameters. Fortunately however, the LARS algorithm enables us to map a dimension to an optimal $\lambda_{n}$. Recall the adaptive LASSO solution path from the LARS in top panel of Figure 1. Transition points (e.g. steps) from 0 to 10 corresponds to a sequence of $\lambda_{n}$’s:

[TABLE]

Note that $\boldsymbol{\hat{\beta}}_{a}(\lambda_{n})=\mathbf{0}$ for $\lambda_{n}>\lambda_{n}(0)$ where $\boldsymbol{\hat{\beta}}_{a}(\lambda_{n})$ is the adaptive LASSO estimator under $\lambda_{n}$. By Theorem 5 in [19],

[TABLE]

where $\hat{df}(\lambda_{n})$ is the number of non-zero elements in $\boldsymbol{\hat{\beta}}_{a}(\lambda_{n})$ and $a_{n}$ is a positive sequence depending on $n$. It is worth mentioning that $\lambda_{n}(m+1)$ is optimum in $[\lambda_{n}(m+1),\lambda_{n}(m))$ by producing the minimum sum of squared errors (SSE) and the smallest model size concurrently.

Also note that the number of steps can exceed the full model size — different steps may have a same model size. Denote $m_{k}$ the last step having a model size $k$, and $m_{k}^{\prime}$ is another step having the same model size. The theorem also showed that

[TABLE]

Theorefore, $\lambda_{n}(m_{k})$ is the overall optimum in $\{\lambda_{n}:\hat{df}(\lambda_{n})=k,\lambda_{n}\in[0,\infty]\}$. So the LARS algorithm enables us to create a one-to-one map between a dimension $k$ and the optimum $\lambda_{n}(m_{k})$,

[TABLE]

It is easy to see that $\lambda_{n}(m_{p_{0}})$ will satisfy condition (A3). Hence, we have the following corollary from Theorem 1.

Corollary 1**.**

Suppose conditions (A1)–(A2) and (A5) hold, then

[TABLE]

where $\mathcal{M}_{r}$ is any $r$-dimensional model.

This result can also be established for adaptive Elastic-Net. Denote $\mathcal{T}_{n}^{*}=\{j:\hat{\beta}_{aej}^{*}\neq 0\}$ an adaptive Elastic-Net estimate of $\mathcal{A}$ using paired bootstrap data.

Corollary 2**.**

Suppose conditions (A1)–(A2) and (A5) hold, then

[TABLE]

where $\mathcal{M}_{r}$ is any $r$-dimensional model.

Proof.

It can be proved by using the techniques for deriving Theorem 1, Corollary 1 and Theorem 2. We bypass here. ∎

We now study the estimation properties for using residual bootstrap data. Denote $\mathcal{T}_{n}^{*}=\{j:\tilde{\beta}_{aej}^{*}\neq 0\}$ an adaptive Elastic-Net estimator of $\mathcal{A}$ using residual bootstrap data.

Theorem 2**.**

Suppose conditions (A1)–(A2) and (A4)–(A5) hold, then

[TABLE]

Moreover, let $\lambda_{n1}^{\prime}$ be another tuning parameter such that the adaptive Elastic-Net estimator under $\lambda_{n1}^{\prime}$ is of dimension $r$, $p_{0}<r<p_{n}$, then

[TABLE]

where $\mathcal{M}_{r}$ is any $r$-dimensional model.

Proofs of Theorem 2 are included in Appendix A. The LARS-EN algorithm for adaptive Elastic-Net estimations is an extension of the LARS algorithm, which shares the same properties of the LARS for deriving Corollaries 1–2. Hence we obtain the following corollary from Theorem 2.

Corollary 3**.**

Suppose conditions (A1)–(A2) and (A5) hold, then

[TABLE]

where $\mathcal{M}_{r}$ is any $r$-dimensional model.

This result can also be established for adaptive LASSO. Denote $\mathcal{A}_{n}^{*}=\{j:\tilde{\beta}_{aj}^{*}\neq 0\}$ an adaptive LASSO estimate of $\mathcal{A}$ using residual bootstrap data.

Corollary 4**.**

Suppose conditions (A1)–(A2) and (A5) hold, then

[TABLE]

where $\mathcal{M}_{r}$ is any $r$-dimensional model.

Proof.

Note that the adaptive LASSO estimator is a special case of the adaptive Elastic-Net estimator with $\lambda_{n2}=0$. Theorem 2 holds automatically for $\mathcal{A}_{n}^{*}$, from which Corollary 4 can be deduced. ∎

Variable selection consistency of the MF procedure can then be deduced from Corollaries 1–4.

Corollary 5**.**

Suppose conditions (A1)–(A2) and (A5) hold. Then the MF procedure is variable selection consistent, e.g.

[TABLE]

where $M_{r^{*}}$ is the model selected from the MF procedure.

Proof.

By definition, $\mathcal{A}_{n}^{*}$ is an adaptive LASSO estimate of $\mathcal{A}$ using paired or residual bootstrap data. It is easy to see that

[TABLE]

Combining with Corollaries 1 or 4,

[TABLE]

Thus the MF procedure for adaptive LASSO can consistently identify the true dimension and true model via selecting the maximum of $\text{MF}_{j},j\in\{1,\dots,p-1\}$, with the highest dimension (if there is a tie). Similarly, Corollary 2 and 3 imply variable selection consistency of the MF procedure for adaptive Elastic-Net. ∎

However, the MF procedure has potential issues in application. In Figure 1 (middle) excluding the full model case, the maximum occurs at dimension 1 instead of 3 although their MFs are both close to 1. In next section, we propose a WMF procedure to tackle this underfitting issue in application.

5 The WMF procedure

5.1 Method and Asymptotic properties

The underfitting issue in MF procedure can be deduced from Corollaries 1–4. Take $\mathcal{A}_{n}^{*}$ for an example. Although it was shown that $\lim_{n\to\infty}P^{*}(\mathcal{A}_{n}^{*}=\mathcal{A}\mid p_{0})=1$, the conditional probability at some underfit dimensions can also reach one, e.g. $\lim_{n\to\infty}P^{*}(\mathcal{A}_{n}^{*}=M_{r}\mid r)=1,0<r<p_{0}$. Note that the tuning parameter leading to an underfit $r$-dimensional estimator, denoted as $\lambda_{n}^{\prime}$, fulfills $\lambda_{n}^{\prime}>\lambda_{n}$. Hence, the convergence rate of $P^{*}(\mathcal{A}_{n}^{*}=M_{r}\mid r)$ at some underfit dimensions can exceed the one at the true dimension. Therefore, the MF procedure would select an underfit model even with a sufficiently large $n$.

In order to fix things, we introduce a weight to the MF procedure. An effective weight should be able to down-weight the underfitting MFs asymptotically, i.e. the weight is able to identify underfit dimensions and its effects does not vanish as $n\to\infty$, without significantly up-weighting the overfitting MFs.

[24] showed that the overall unconditional (on $\mathbf{y}$) expected squared prediction error for the OLS estimator of $\boldsymbol{\beta}_{0}$ under model $\alpha$ is

[TABLE]

where $p_{\alpha}$ indicates the size of $\alpha$, $\Delta_{\alpha,n}=\boldsymbol{\beta}_{0}^{T}\mathbf{X}^{T}(\mathbf{I}-\mathbf{P}_{\alpha})\mathbf{X}\boldsymbol{\beta}_{0}/n$,

$\mathbf{P}_{\alpha}=\mathbf{X}_{\alpha}(\mathbf{X}_{\alpha}^{T}\mathbf{X}_{\alpha})^{-1}\mathbf{X}_{\alpha}^{T}$, $\mathbf{X}_{\alpha}$ is a sub-matrix of $\mathbf{X}$ whose columns are indexed by the components of $\alpha$ and $\mathbf{I}$ is an identity matrix.

When $\alpha$ is a true or overfit model, it has $\mathbf{X}\boldsymbol{\beta}_{0}=\mathbf{X}_{\alpha}\boldsymbol{\beta}_{\alpha}$ and thus

[TABLE]

However, if $\alpha$ is an underfit model, then $\Delta_{\alpha,n}>0$ for any fixed $n$. He further assumed that

[TABLE]

which is argued in the paper to be a minimal type of asymptotic model identifiability condition. Under assumption (5.3) and by (5.1)–(5.2),

[TABLE]

where $\nu$ is an underfit model and $\kappa$ is a true or overfit model. By (5.4) a formula inversely proportional to $T_{\alpha,n}$ will be an ideal choice for the weight.

[25] proposed such a formula for estimating the posterior probability of the model size given the data

[TABLE]

where $\hat{T}_{n}(j)$ is an estimate of $T_{\alpha,n}$ using a $j$-dimensional model and $c,1\leq c\leq 2$, is a constant. We use the multi-fold CV for $\hat{T}_{n}(j)$ and define

[TABLE]

Figure 1 (bottom) shows the effect of weights in Example 1, which heavily punish underfitting MFs and have little effect on true and overfitting MFs. The WMF procedure then selects the dimension $r^{*}$ and model $M_{r^{*}}$ s.t.

[TABLE]

Recall that $\text{MF}_{j}/B$ is a bootstrap version estimate of the posterior probability of model $M_{j}$ given the data and dimension, i.e. $P(M_{j}\mid\mathbf{y},j)$, along with (5.6) it has

[TABLE]

Note that BIC is a Laplace approximation to $P(M_{j}\mid\mathbf{y})$ under a flat prior assumption and is variable selection consistent for adaptive LASSO [11, 12], but no convergence rate has been studied. Simulation studies in Section 6 show that BIC has a much slower empirical convergence rate than the WMF procedure.

Next we show properties of the multi-fold CV using adaptive LASSO or adaptive Elastic-Net estimators. Then variable selection consistency of the WMF procedure can be established. Let $K$ be a fixed integer and suppose $n=Kt$. In multi-fold CV, one randomly divides a sample of $n$ observations into $K$ mutually exclusive subgroups $s_{1},\dots,s_{K}$ with each subgroup containing $t$ observations, and selects the model by minimizing the following sum of squared errors

[TABLE]

where $\boldsymbol{\hat{\beta}}_{s_{i}^{c},\mathcal{M}}$ is an adaptive LASSO or adaptive Elastic-Net estimator under model $\mathcal{M}$ using samples not in $s_{i}$. Let $\alpha$ and $\alpha^{\prime}$ be the true or overfit models and $\nu$ be an underfit model. We assume following condition for asymptotic studies of the multi-fold CV procedure.

(A6) $\sup_{t\to\infty}\sup_{s_{i}}\|t^{-1}\mathbf{X}_{s_{i},\mathcal{M}}^{T}\mathbf{X}_{s_{i},\mathcal{M}}-\mathbf{V}_{\mathcal{M}}\|=o(1),$ where $\mathbf{V}_{\mathcal{M}}$ is a positive definite matrix.

Theorem 3**.**

*Suppose conditions (A1)–(A2) and (A5)–(A6) hold, then

1. the multi-fold CV for adaptive LASSO or adaptive Elastic-Net satisfies*

[TABLE]

2. model $M_{r^{*}}$ selected from the WMF procedure fulfills

[TABLE]

Proofs of Theorem 3 are included in Appendix A. Denote $r^{\prime}$ an underfit dimension. The ratio of ${\text{WMF}_{p_{0}}\over\text{WMF}_{r^{\prime}}}$ is exponentially proportional to the bias term, ${d\over 2c\sigma^{2}}\|\boldsymbol{\beta}_{0\mathcal{M}_{r^{\prime}}^{c}}\|^{2}$, which is larger than 0 and does not fade as $n\to\infty$. This guarantees a good finite sample performance of the WMF procedure and a fast vanishing rate of its underfitting issues, which will be confirmed in simulation studies in Section 6.

5.2 Computation

In adaptive Elastic-Net, $\lambda_{n2}$ takes the same value in Elastic-Net for calculating the weights $\omega_{j}$’s, where the tuning parameters are chosen by minimizing the two-dimensional BIC [7]. Then computational efforts remain the same for adaptive LASSO and adaptive Elastic-Net, which are to compute a full solution path against $\lambda_{n}$’s or $\lambda_{n1}^{+}$’s. Computational complexity of creating an entire adaptive LASSO solution path is of order $O(np_{n}^{2})$ [6]. It is of order $O(np_{n}^{2}+p_{n}^{3})$ for adaptive Elastic-Net[7]. Since the optimal value often occurs at an early stage, we could stop the algorithms after $m,m<p_{n},$ steps. In this case, the computational cost reduces to $O(nm^{2})$ for adaptive LASSO and $O(m^{3}+nm^{2})$ for adaptive Elastic-net.

Computational cost of a WMF procedure is then $B$ times the cost of computing an adaptive LASSO or adaptive Elastic-Net solution path.

6 Empirical studies

We now investigate empirical performances of the WMF procedure and show it outperforms the BIC, EBIC, GIC, SS, Cp, and 1se-CV (which is often recommended for variable selection) in a wide range of situations for both adaptive LASSO and adaptive Elastic-Net. The Cp did very poor in all scenarios, thus is excluded in the presentation.

In all simulations, data were generated from

[TABLE]

where $\mathbf{x}_{i}\overset{iid}{\sim}N_{p_{n}}(\mathbf{0},\boldsymbol{\Sigma})$ and $\varepsilon_{i}\overset{iid}{\sim}N(0,1)$. Let $p_{n}=O(n^{\kappa})$ for some constant $\kappa$, $0\leq\kappa<1$, $n=100,300,500$. Results were averaged over 100 times of replications.

6.1 Simulations of the adaptive LASSO WMF procedure

Three scenarios were designed for the adaptive LASSO WMF procedure. In each scenario, $\boldsymbol{\Sigma}(i,j)=0.3^{|i-j|}$ and $\sigma=3$.

Scenario 1: Fixed low dimension and moderate proportion of true covariates. More specifically, set $p_{n}=10$ and $\boldsymbol{\beta}=(3,1.5,0,0,2,0,\dots)_{10}^{T}$. Then the proportion of true covariates is 0.3, and the signal to noise ratios (SNR) are respectively 2.03, 2 and 1.98 for various $n$.

Scenario 2: Low dimension, moderate proportion of true covariates and weak signals for some true covariates. Specifically, set $p_{n}=O(\sqrt{n})$, then $p_{n}$ equals to 10, 17, 22 accordingly. Let $p_{0}$ grow with $n$ as follows. Initially $p_{0}=3$ and $\boldsymbol{\beta}=(3,1.5,0,0,2,0,\dots)^{T}$. Afterwards, $p_{0}$ increases by 1 for every 40-unit increment in $n$ and the new element equals to 1. As a result, the proportions of true covariates are respectively 0.3, 0.47, and 0.59, and the SNRs are 2, 2.85 and 3.69.

Scenario 3: High dimension, sparse proportion of true covariates and relatively large signals for all true covariates. In detail, set $p_{n}=O(n^{3/4})$, then $p_{n}$ equals to 32, 72, 106 accordingly. Let $p_{0}$ grow in the same manner as in scenario 2, but the new elements equal to 2. Accordingly, the proportions of true covariates are 0.09, 0.11 and 0.12, and the SNRs are 2, 5.07, and 8.5.

Paired bootstrapping was used in the adaptive LASSO WMF procedure. Simulation results are summarized in Figures 2–4. In all scenarios, the proposed method has the highest degree of accuracy in identifying the true model and also enjoys a much faster convergence rate than other compared methods. The WMF procedure has an underfitting issue which vanishes quickly as $n$ increases. Other methods (except for the SS) however suffer from an overfitting issue. The sparser the model is, the more serious the issue tends to be. Performance of the SS relies on particular specifications of several unknown parameters. Although we have followed instructions in [17] for setting those parameters throughout the simulations, its performance remains erratic and unsatisfactory.

Simulations for using residual bootstrapping in the adaptive LASSO WMF procedure were also conducted. The results are presented in Appendix B, which are similar to those in above paired bootstrapping simulations.

6.2 Simulations of the adaptive Elastic-Net WMF procedure

We also designed three scenarios for the adaptive Elastic-Net WMF procedure, each of which mimics a typical structure in applications. Since the adaptive Elastic-Net fits data with grouping effects, in following simulations true covariates will be added in blocks with size 3. The SS is excluded due to its poor performance.

Scenario 4: Low dimension, moderate proportion of true covariates, weak signals for some true covariates and moderate correlations between covariates. More specifically, let $\boldsymbol{\Sigma}(i,j)=0.5^{|i-j|}$, $\sigma=3$, and $p_{n}=O(\sqrt{n})$. Initially we have one block of true covariates, then $p_{0}=3$. Elements of $\boldsymbol{\beta}$ in the block equal to 2, the rest are 0. Afterwards, we add 1 block of true covariates for every 200-unit increment in $n$ and the new elements equal to 1. Respectively, the proportions of true covariates are 0.3, 0.35 and 0.41, and the SNRs are 2.45, 3.72 and 3.67.

Scenario 5: High dimension, sparse proportion of true covariates, relatively large signals for all true covariates and moderate correlations between covariates. In detail, let $\boldsymbol{\Sigma}(i,j)=0.5^{|i-j|}$, $\sigma=5$, and $p_{n}=O(n^{3/4})$. Initially set $p_{0}=6$. Then true covariates follow the same adding scheme as in scenario 4. All non-zero elements in $\boldsymbol{\beta}$ equal to 2. Respectively, the proportions of true covariates are 0.19, 0.13 and 0.11, and the SNRs are 1.79, 3.09 and 3.51.

Scenario 6: High dimension, sparse proportion of true covariates, relatively large signals for all true covariates and high correlations between grouped covariates. Specifically, let $\sigma=5$ and $p_{n}=O(n^{3/4})$. True covariates follow the same adding scheme as in Scenario 5, all non-zero elements in $\boldsymbol{\beta}$ equal to 2. Moreover, true covariates within each block have correlations almost 1, while true covariates between the blocks have correlation 0. All noise covariates are i.i.d from $N(0,1)$. Respectively, the proportions of true covariates are 0.19, 0.13 and 0.11, and the SNRs are 2.84, 4.35 and 5.75.

Residual bootstrapping was used in the adaptive Elastic-Net WMF procedure. Simulation results are summarized in Figures 5–7. In scenarios 4 and 5, the proposed method has the best performance over other compared methods: on average the highest degree of accuracy in indentifying the true model; a faster convergence rate; the underfitting issue vanishes quickly. On the other hand, other methods suffer from an overfitting issue. The sparser the model is, the more serious the issue tends to be. In scenario 6, all methods do equally well because the adaptive Elastic-Net well fit the data with highly grouped effects.

Simulation results for using paired bootstrapping in the adaptive Elastic-Net WMF procedure are presented in Appendix B, which are similar to those in above residual bootstrapping simulations.

6.3 Classification analysis of the leukaemia data

We now demonstrate the WMF procedure in a real data application. The leukaemia data [26] contains $p_{n}=7129$ genes and $n=72$ samples. We have 38 out of the 72 samples from the training dataset with 27 ALL’s (acute lymphoblastic leukaemia) and 11 AML’s (acute myeloid leukaemia). The remaining 34 samples are from the test dataset with 20 ALL’s and 14 AML’s. The goal of this analysis is to identify a subset of genes that can accurately predict the type of leukaemia for future data. Similar to [7], we coded the type of leukaemia as a binary response variable, denoted as $y$, and defined the classification function as $I(\hat{y}>0.5)$, where $I(\cdot)$ is the indicator function.

To improve computational efficiency, we selected 1000 candidate genes as the predictors using the sure independence screening (SIS) procedure [27]. The adaptive LASSO and adaptive Elastic-Net were then applied to explore the data. The screening and variable selection were carried out on the training dataset, while classification errors were examined on the test dataset. Both the LARS and LARS-EN algorithms were stopped after 200 steps of estimation to further reduce the computational costs. Note that since the optimal steps selected by various types of methods are much smaller than the stopping step, this strategy will not affect the variable selection.

Classification results are summarized in Tables 1–2. For adaptive LASSO, although the Cp, CV and BIC have obtained the minimal classification errors for both training and test datasets, the WMF has classification errors close to the minimum using the least number of genes. For adaptive Elastic-Net, the WMF has the minimal classification errors for both training and test datasets using the least number of genes. Thus we conclude that the WMF procedure is able to find the set of “important” genes that can largely improve the prediction accuracy.

7 Extensions

Here we investigate extensions of the WMF procedure to GLMs, which has the following generic density fuction [28]

[TABLE]

[6] had extended the adaptive LASSO to GLMs. Its estimator, $\boldsymbol{\hat{\beta}}_{a}$, is obtained by maximizing the penalized log-likelihood,

[TABLE]

where $\hat{w}_{j}=1/|\tilde{\beta}_{j}|^{\gamma}$, $\gamma>0$ and $\boldsymbol{\tilde{\beta}}=(\tilde{\beta}_{1},\dots,\tilde{\beta}_{p})^{T}$ is the maximum likelihood estimator. Under certain regularity conditions, $\boldsymbol{\hat{\beta}}_{a}$ was shown to enjoy the oracle properties .

The generalization of Multi-fold CV to GLMs is straightforward [24]. Define,

[TABLE]

where $Q(\cdot,\cdot)$ is a loss function, $\mathbf{\hat{y}}_{s_{i}^{c},\alpha}$ is the prediction of $\mathbf{y}_{s_{i}}$ under model $\alpha$ using samples not in $s_{i}$.

Then we can extend the WMF procedure to GLMs for adaptive LASSO. In this case, we draw $B$ paired bootstrap samples in step 1 of Algorithm 1. Note that the LARS algorithm does not fit for GLMs, but we can use the coordinate descent algorithm [29] instead, which generates a solution path similar to the LARS. Hence in step 2, we use the coordinate descent algorithm to fit each bootstrap data. The rest remain the same. Asymptotic properties of the adaptive LASSO WMF procedure for GLMs can also be established by using some similar techniques for showing Theorem 1 in this paper and Theorem 4 in [6].

We demonstrate this extension through one simple example, where binary responses were generated from the logistic regression model

[TABLE]

where $\mathbf{x}_{i}\overset{iid}{\sim}N_{10}(\mathbf{0},\boldsymbol{\Sigma})$, $\boldsymbol{\Sigma}(i,j)=0.3^{|i-j|}$, and $\boldsymbol{\beta}=(3,1.5,0,0,2,0,\dots)_{10}^{T}$. Simulation results were averaged over 100 times of replications and summarized in Figure 8. It shows that the WMF procedure is much more accurate in variable selection and also enjoys a faster convergence rate than other compared methods.

Extension of the adaptive Elastic-Net WMF procedure to GLMs is similar. Define the adaptive Elastic-Net estimator for GLMs as

[TABLE]

where $w_{j}=|\hat{\beta}_{ej}|^{-\gamma}$, $\gamma>0$ and $\boldsymbol{\hat{\beta}}_{e}=(\hat{\beta}_{e1},\dots,\hat{\beta}_{ep_{n}})^{T}$ is defined in (7.1) with $\hat{w}_{j}=1$ for all $j$’s. The rest follow the same procedures for extension of the adaptive LASSO WMF procedure.

8 Ultra-high dimensional data

In this section, we discuss applications of the WMF procedure to ultra-high dimensional data in which $p_{n}>n$. [27] proposed the sure independence screening (SIS) method for ultra-high dimensional data to reduce their dimensionality to a moderate scale, $d_{n}$, s.t. $d_{n}<n$. Afterwards a lower dimensional estimation method such as the SCAD can be applied to the reduced data. This process is called SIS+SCAD. Under some regularity conditions, they showed that the SIS has an exponentially small probability to omit true features and the SIS+SCAD retains the oracle properties if $d_{n}=o_{p}(n^{1/3})$. By replacing the SCAD with adaptive Elastic-Net, the new procedure is refered to as SIS+AEnet [8], which holds the oracle properties if $d_{n}=O_{p}(n^{\varrho}),0\leq\varrho<1$. Here we recommend to combine SIS with the WMF procedure when $p_{n}>n$. We first use the SIS to reduce the dimensionality to $d_{n},d_{n}<n$, and then apply the WMF procedure to the reduced data. We call this procedure SIS+WMF.

Corollary 6**.**

Suppose conditions for Theorem 1 in [27] and Theorem 3 in this paper hold. Let $d_{n}=n^{\varrho},0\leq\varrho<1$. Then the SIS+WMF procedure is variable selection consistent.

Note that Corollary 6 is a direct conclusion of Theorem 1 in [27] and Theorem 3 in this paper.

9 Discussion

We proposed a prediction-weighted maximal frequency procedure to estimate the amount of regularization for adaptive LASSO and adaptive Elastic-Net. Asymptotic properties were studied with a diverging $p_{n}$.

Central idea of the WMF procedure is the importance of conditioning on dimension, which mitigates overfitting. Underfitting can then be handled by using prediction-based weights estimated by multi-fold cross-validation. This simple recipe can also be applied to other regularization methods, say the SCAD and fused LASSO, making the WMF procedure a unified model selection criterion in regularization problems. However, asymptotic properties have yet to be studied, which will be a future topic.

Appendix A Proofs

Proof of Lemma 1.

Assume $\mid\beta_{i}\mid>\mid\beta_{j}\mid$ and $\mid\beta_{i}\mid-\mid\beta_{j}\mid=m\sigma,\ m>0$. We have 4 cases for $\beta_{i},\beta_{j}$

[TABLE]

Let $Z_{i}=\beta_{i}+{x}_{i}^{T}\varepsilon\sim N(\beta_{i},\sigma^{2})$ and $Z_{j}=\beta_{j}+{x}_{j}^{T}\varepsilon\sim N(\beta_{j},\sigma^{2})$. We have

[TABLE]

Consider case 1: $\beta_{j}\geq 0$ and $\beta_{i}=\beta_{j}+m\sigma$, $m>0$.

Let $k$ be a positive constant. The point $\beta_{j}+k\sigma$ separates the domain of $Z_{i}$ and $Z_{j}$ into two parts: $(-\infty,\;\beta_{j}+k\sigma]$ and $(\beta_{j}+k\sigma,\;\infty]$. The cumulative probabilities of $Z_{i}$ and $Z_{j}$ in first part of the domain are respectively

[TABLE]

The probability $P(\mid Z_{i}\mid>\mid Z_{j}\mid)$ can then be calculated from

[TABLE]

After some simple deductions, we get,

[TABLE]

If $m=0$ i.e. $\mid\beta_{i}\mid=\mid\beta_{j}\mid$, from (A.1) we have

[TABLE]

However if $m>0$ i.e. $\mid\beta_{i}\mid>\mid\beta_{j}\mid$,

[TABLE]

Since $m,k,\beta_{j},\sigma>0$, we have

[TABLE]

Note that two integrals in (A.2) have equal length of the integral intervals. Moreover the integral function is an monotonically decreasing function of $x$ for $x\geq 0$, and monotonically increasing for $x<0$. Hence

[TABLE]

Combining (A.2) with (A.3), we get

[TABLE]

Other three cases can be proved in the same way. We avoid the repetitions here. ∎

Proof of Theorem 1.

By [8], $\boldsymbol{\hat{\beta}}_{a}$ enjoys the oracle properties under certain regularity conditions. And $\boldsymbol{\hat{\beta}}_{a}^{*}$ is a paired bootstrap analog of $\boldsymbol{\hat{\beta}}_{a}$ by replacing $(\mathbf{X},\mathbf{y})$ with $(\mathbf{X}^{*},\mathbf{y}^{*})$ in estimation. To simplify notations in the proof, we drop the subscript ‘$a$’ in $\boldsymbol{\hat{\beta}}_{a}$ and $\boldsymbol{\hat{\beta}}_{a}^{*}$.

By the KKT regularity conditions, $\boldsymbol{\hat{\beta}}^{*}$ is the unique solution of adaptive LASSO given $(\mathbf{X}^{*},\mathbf{y}^{*})$ if

[TABLE]

where $\mathbf{X}_{j}^{*}$ is the $j$th column of $\mathbf{X}^{*}$ and

[TABLE]

Let $\mathbf{\tilde{s}}_{\mathcal{A}}=(\omega_{j}sgn(\hat{\beta}_{j}),j\in\mathcal{A})^{T}$ and $\boldsymbol{\hat{\beta}}_{\mathcal{A}}^{*}=(\mathbf{X}_{\mathcal{A}}^{*T}\mathbf{X}_{\mathcal{A}}^{*})^{-1}(\mathbf{X}_{\mathcal{A}}^{*T}\mathbf{y}^{*}-\lambda_{n}\mathbf{\tilde{s}}_{\mathcal{A}})$. We show that $(\boldsymbol{\hat{\beta}}_{\mathcal{A}}^{*},\mathbf{0})$ satisfies (A.4) with probability tending to 1, which is equivalent to prove

[TABLE]

where the first inequation implies $sgn(\boldsymbol{\hat{\beta}}_{\mathcal{A}}^{*})=sgn(\boldsymbol{\hat{\beta}}_{\mathcal{A}})$.

Note that $\omega_{j}=|\tilde{\beta}_{j}|^{-\gamma}$, where $\boldsymbol{\tilde{\beta}}=(\tilde{\beta}_{1},\dots,\tilde{\beta}_{p_{n}})^{T}$ is an OLS or best ridge estimate of $\boldsymbol{\beta}_{0}$,

[TABLE]

By Theorem 3.1 in [8],

[TABLE]

under assumption that $\lim_{n\to\infty}{\lambda_{n2}\over\sqrt{n}}=0$. It is satisfied automatically for the OLS estimate.

Denote $\mathbf{x}_{i\mathcal{A}}^{*}$ the $i$th row of $\mathbf{X}_{\mathcal{A}}^{*}$, and $\otimes$ the element-wise product. We have

[TABLE]

Hence under conditions (A1) and (A5),

[TABLE]

Let $\psi=\min_{j\in\mathcal{A}}|\beta_{0j}|$, $\tilde{\psi}=\min_{j\in\mathcal{A}}|\tilde{\beta}_{j}|$ and $\hat{\psi}=\min_{j\in\mathcal{A}}|\hat{\beta}_{j}|$. Under conditions (A1)–(A3) and (A5), the first inequation in (A.5) can be proved by

[TABLE]

where

[TABLE]

By (A.6), it has

[TABLE]

Similarly,

[TABLE]

By Theorem 3.1 in [8],

[TABLE]

For proof of the second inequation in (A.5), it suffices to show

[TABLE]

Since

[TABLE]

it follows that

[TABLE]

where $\kappa$, $0<\kappa<1$, is a constant.

For $B_{1}$,

[TABLE]

where $p_{\mathcal{A}^{c}}$ indicates the size of $\mathcal{A}^{c}$. By (A.6), $\forall j\in\mathcal{A}^{c}$, $E|\tilde{\beta}_{j}|^{2}\leq E\|\boldsymbol{\tilde{\beta}}-\boldsymbol{\beta}_{0}\|^{2}=O_{p}\left({p_{n}\over n}\right)$, which indicates $|\tilde{\beta}_{j}|\leq O_{p}\left({p_{n}\over n}\right)^{1/2}$. Then under condition (A3), $B_{1}$ fulfills

[TABLE]

Also since

[TABLE]

we have for $B_{2}$,

[TABLE]

Hence (A.5) is proved. We have shown that $\boldsymbol{\hat{\beta}}^{*}=(\boldsymbol{\hat{\beta}}_{\mathcal{A}}^{*},\mathbf{0})$ and $sgn(\boldsymbol{\hat{\beta}}_{\mathcal{A}}^{*})=sgn(\boldsymbol{\hat{\beta}}_{\mathcal{A}})$ with probability tending to 1, where $\boldsymbol{\hat{\beta}}^{*}$ is the adaptive LASSO estimate using paired bootstrap data. Also it can be deduced from (A.7) that

$P(\min_{j\in\mathcal{A}}|\hat{\beta}_{j}|>0)\to 1$. To sum up, we get $\lim_{n\to\infty}P^{*}(\mathcal{A}_{n}^{*}=\mathcal{A}\mid\lambda_{n})=1$.

We now prove $\lim_{n\to\infty}P^{*}(\mathcal{A}_{n}^{*}=\mathcal{M}_{r}\mid\lambda_{n}^{\prime})<1$, where $\mathcal{M}_{r}$ is any $r$-dimensional model, $p_{0}<r<p_{n}$, and $\lambda_{n}^{\prime}$ is a tuning parameter such that the adaptive LASSO estimator under $\lambda_{n}^{\prime}$ is of dimension $r$. Then $\lambda_{n}^{\prime}<\lambda_{n}$, hence $\lambda_{n}^{\prime}/{\sqrt{n}}\to 0$. If it also satisfies $\lim_{n\to\infty}{\lambda_{n}^{\prime 2}n^{(1-\varrho)(1+\gamma)-1}\over n}\to\infty$, we would have $P^{*}(\mathcal{A}_{n}^{*}=\mathcal{A}\mid\lambda_{n}^{\prime})=1$ based on previous proof, which contradicts with the definition of $\lambda_{n}^{\prime}$. Therefore,

[TABLE]

To prove $\lim_{n\to\infty}P^{*}(\mathcal{A}_{n}^{*}=\mathcal{M}_{r}\mid\lambda_{n}^{\prime})<1$, by the KKT regularity conditions it suffices to show

[TABLE]

or equivalently

[TABLE]

Following previous proof, we get

[TABLE]

However,

[TABLE]

as $n\to\infty$. Similarly, $\lim_{n\to\infty}B_{2}\not\to 0$. Then (A.8) holds. ∎

Lemma 2**.**

Suppose conditions (A1) and (A5) hold and $\lim_{n\to\infty}\lambda_{n2}/\sqrt{n}=0$ in ridge estimates. Then,

[TABLE]

Proof.

Assume $\boldsymbol{\hat{\beta}}$ is a ridge estimate of $\boldsymbol{\beta}_{0}$,

[TABLE]

By (A.6), E\|\boldsymbol{\hat{\beta}}-\boldsymbol{\beta}_{0}\|^{2}\leq O_{p}\big{(}{p_{n}\over n}\big{)}. Calculate centered residuals $\boldsymbol{\hat{\varepsilon}}$,

[TABLE]

where each entry of $\boldsymbol{\bar{\varepsilon}}_{0}$, marked as $\bar{\varepsilon}_{0}$, is the mean of $\boldsymbol{\hat{\varepsilon}}_{0}$. Denote $\boldsymbol{\varepsilon}^{*}=(\varepsilon_{1}^{*},\dots,\varepsilon_{n}^{*})^{T}$ an i.i.d bootstrap sample from the empirical distribution that puts mass $n^{-1}$ on each entry of $\boldsymbol{\hat{\varepsilon}}$.

By definition, we have

[TABLE]

[TABLE]

and

[TABLE]

In above equation,

[TABLE]

Moreover, by the sum of squares inequality,

[TABLE]

Hence,

[TABLE]

Let

[TABLE]

where $\bar{\varepsilon}={1\over n}\sum_{i=1}^{n}\varepsilon_{i}$. We now prove $s_{n}\to\sigma_{n}$ asymptotically.

Note that

[TABLE]

with probability 1.

And by the sum of squares inequality,

[TABLE]

Then $\lim_{n\to\infty}s_{n}^{2}=\sigma^{2}$ with probability 1. ∎

Proof of Theorem 2.

Let $(\mathbf{X},\mathbf{y}^{*})$ be a residual bootstrap sample, where $\mathbf{y}^{*}=\mathbf{X}\boldsymbol{\hat{\beta}}+\boldsymbol{\varepsilon}^{*}$ and $\boldsymbol{\hat{\beta}}$ is the ridge estimator. Define

[TABLE]

where we dropped the subscript ‘ae’ in $\boldsymbol{\tilde{\beta}}_{ae}^{*}$ for simplicity.

Let

[TABLE]

we prove $((1+{\lambda_{n2}\over n})\boldsymbol{\tilde{\beta}}_{\mathcal{A}}^{*},\mathbf{0})$ is the solution to (A.9) with probability tending to 1. By the KKT regularity conditions, this suffices to show

[TABLE]

or equivalently

[TABLE]

Note that $\omega_{j}=|\hat{\beta}_{ej}|^{-\gamma}$ where $\boldsymbol{\hat{\beta}}_{e}=(\hat{\beta}_{e1},\dots,\hat{\beta}_{ep_{n}})^{T}$ is the Elastic-Net estimator defined in (2.2). By Theorem 3.1 in [8],

[TABLE]

under condition (A4).

Let $\psi=\min_{j\in\mathcal{A}}|\beta_{0j}|$ and $\tilde{\psi}=\min_{j\in\mathcal{A}}|\hat{\beta}_{ej}|$. Then

[TABLE]

By (A.11) under condition (A4),

[TABLE]

Also by (A.11) $\forall j\in\mathcal{A}^{c}$, $E|\hat{\beta}_{ej}|^{2}\leq E\|\boldsymbol{\hat{\beta}}_{e}-\boldsymbol{\beta}_{0}\|^{2}=O_{p}\left({p_{n}\over n}\right)$, which indicates $|\hat{\beta}_{ej}|\leq O_{p}\left({p_{n}\over n}\right)^{1/2}$. Hence

[TABLE]

Note that

[TABLE]

By (A.6),

[TABLE]

We now study $E^{*}\|\boldsymbol{\hat{\beta}}_{\mathcal{A}}-\boldsymbol{\tilde{\beta}}_{\mathcal{A}}^{*}\|^{2}$. Let

[TABLE]

By using the same arguments for deriving (6.3) in [8], we can easily show

[TABLE]

On the other hand,

[TABLE]

by Lemma 2,

[TABLE]

By assembling (A.12)–(A.14), we get

[TABLE]

And

[TABLE]

Then under conditions (A1)–(A2) and (A4)–(A5),

[TABLE]

Hence (A.10) is proved. So far we have shown that $\boldsymbol{\tilde{\beta}}^{*}=((1+{\lambda_{n2}\over n})\boldsymbol{\tilde{\beta}}_{\mathcal{A}}^{*},\mathbf{0})$ with probability tending to 1, where $\boldsymbol{\tilde{\beta}}^{*}$ is the adaptive Elastic-Net estimate using residual bootstrap data. To prove $\lim_{n\to\infty}P^{*}(\mathcal{T}_{n}^{*}=\mathcal{A}\mid\lambda_{n1}^{+})=1$, we still need to show that $P(\min_{j\in\mathcal{A}}|\tilde{\beta}_{j}^{*}|>0)\to 1$.

Let $\hat{\psi}=\min_{j\in\mathcal{A}}|\hat{\beta}_{j}|$ and $\tilde{\psi}^{*}=\min_{j\in\mathcal{A}}|\tilde{\beta}_{j}^{*}|$. By (A.6),

[TABLE]

Hence $P(\hat{\psi}>\psi/2)\to 1$ as $n\to\infty$ where $\psi>0$. Under condition (A4),

[TABLE]

which indicates $P(\tilde{\psi}^{*}>\hat{\psi}/2)\to 1$ as $n\to\infty$. To sum up, $\lim_{n\to\infty}P(\tilde{\psi}^{*}>\psi/4)=1$. Thus $\lim_{n\to\infty}P^{*}(\mathcal{T}_{n}^{*}=\mathcal{A}\mid\lambda_{n1}^{+})=1$ is proved.

We now prove $\lim_{n\to\infty}P^{*}(\mathcal{T}_{n}^{*}=\mathcal{M}_{r}\mid\lambda_{n1}^{\prime})<1$, where $\mathcal{M}_{r}$ is any $r$-dimensional model, $p_{0}<r<p_{n}$, and $\lambda_{n1}^{\prime}$ is a tuning parameter such that the adaptive Elastic-Net estimator under $\lambda_{n1}^{\prime}$ is of dimension $r$. Then $\lambda_{n1}^{\prime}<\lambda_{n1}^{+}$, hence $\lambda_{n1}^{\prime}/\sqrt{n}\to 0$. If it also satisfies $\lim_{n\to\infty}{\lambda_{n1}^{\prime 2}n^{(1-\varrho)(1+\gamma)-1}\over n}\to\infty$, we would have $P^{*}(\mathcal{T}_{n}^{*}=\mathcal{A}\mid\lambda_{n1}^{\prime})=1$ based on previous proof, which contradicts with the definition of $\lambda_{n1}^{\prime}$. Therefore,

[TABLE]

To prove $\lim_{n\to\infty}P^{*}(\mathcal{T}_{n}^{*}=\mathcal{M}_{r}\mid\lambda_{n1}^{\prime})<1$, by the KKT regularity conditions it suffices to show

[TABLE]

or equivalently

[TABLE]

By following the same arguments for showing (A.10), we get

[TABLE]

∎

Lemma 3**.**

Suppose conditions (A1), (A5) and (A6) hold. Denote $\alpha$ an overfit model including the true model, the adaptive Elastic-Net estimate $\boldsymbol{\hat{\beta}}_{s_{i}^{c},\alpha}$ from the multi-fold CV then satisfies

[TABLE]

where the adaptive LASSO estimate is a special case with $\lambda_{n2}=0$.

Proof.

Here we provide a proof for the adaptive LASSO estimator. The adaptive Elastic-Net estimator can be proved by using the same arguments for deriving Theorem 3.1 in [8] and the strategies in below.

The adaptive LASSO estimator from the multi-fold CV is

[TABLE]

which satisfies

[TABLE]

Hence,

[TABLE]

The last equation holds because $\lambda_{n1}^{\prime}$ continuously decreases from $\lambda_{n1}^{+}$ to 0 as $\alpha$ changes from the true model to full model. ∎

Proof of Theorem 3.

We integrate the proof for adaptive Elastic-Net and adaptive LASSO. Denote $\alpha$ an overfit model including the true model. The $\text{MCV}_{\alpha}$ is

[TABLE]

By Lemma 3, the second term in (A.15) satisfies

[TABLE]

Hence,

[TABLE]

The third term in (A.15) fulfills

[TABLE]

Hence,

[TABLE]

By substituting (A.16)–(A.17) to (A.15), we obtain

[TABLE]

Let $\alpha$ and $\alpha^{\prime}$ be two overfit models including the true model, then

[TABLE]

We now consider an underfit model $\nu$. The $\text{MCV}_{\nu}$ is

[TABLE]

Let $\boldsymbol{\hat{\beta}}_{\nu}$ be an adaptive LASSO or adaptive Elastic-Net estimator under $\nu$. The second term in (A.18) satisfies

[TABLE]

For the third term in (A.18),

[TABLE]

Hence,

[TABLE]

By substituting (A.19)–(A.20) to (A.18), we get

[TABLE]

If $\alpha$ is an overfit model and $\nu$ is an underfit model, we have

[TABLE]

So the first part is proved. We then combine it with Corollaries 1–4. For any $r$, $p_{0}<r<p_{n}$,

[TABLE]

And for any $r^{\prime}$, $0<r^{\prime}<p_{0}$,

[TABLE]

Then model selection consistency of the WMF procedure can be deduced from (A.22)–(A.23). ∎

Appendix B Additional simulation results