Open AccessArticle

The One Standard Error Rule for Model Selection: Does It Work?

Yuchen Chen

¹ and

Yuhong Yang

^2,*

Carlson School of Management, University of Minnesota, Minneapolis, MN 55455, USA

School of Statistics, University of Minnesota, Minneapolis, MN 55455, USA

Author to whom correspondence should be addressed.

Stats 2021, 4(4), 868-892; https://doi.org/10.3390/stats4040051

Submission received: 25 September 2021 / Revised: 28 October 2021 / Accepted: 1 November 2021 / Published: 5 November 2021

(This article belongs to the Special Issue Re-sampling Methods for Statistical Inference of the 2020s)

Download

Browse Figures

Versions Notes

Abstract

Previous research provided a lot of discussion on the selection of regularization parameters when it comes to the application of regularization methods for high-dimensional regression. The popular “One Standard Error Rule” (1se rule) used with cross validation (CV) is to select the most parsimonious model whose prediction error is not much worse than the minimum CV error. This paper examines the validity of the 1se rule from a theoretical angle and also studies its estimation accuracy and performances in applications of regression estimation and variable selection, particularly for Lasso in a regression framework. Our theoretical result shows that when a regression procedure produces the regression estimator converging relatively fast to the true regression function, the standard error estimation formula in the 1se rule is justified asymptotically. The numerical results show the following: 1. the 1se rule in general does not necessarily provide a good estimation for the intended standard deviation of the cross validation error. The estimation bias can be 50–100% upwards or downwards in various situations; 2. the results tend to support that 1se rule usually outperforms the regular CV in sparse variable selection and alleviates the over-selection tendency of Lasso; 3. in regression estimation or prediction, the 1se rule often performs worse. In addition, comparisons are made over two real data sets: Boston Housing Prices (large sample size n, small/moderate number of variables p) and Bardet–Biedl data (large p, small n). Data guided simulations are done to provide insight on the relative performances of the 1se rule and the regular CV.

Keywords:

subsampling; tuning parameter selection; estimation accuracy; regression estimation; variable selection

1. Background

Resampling and subsampling methods were widely utilized in many statistical and machine learning applications. In high-dimensional regression learning, for instance, bootstrap and subsampling play important roles in model selection uncertainty quantification and other model selection diagnostics (see, e.g., [1] for references). Cross-validation (CV) (as a subsampling tool) and closely related methods were proposed to assess quality of variable selection performance in terms of F- and G-measures ([2]) and to determine variable importances ([3]).

In this paper, we focus on the examination of cross-validation as used for model selection. A core issue is how to properly quantify variabilities in subsamplings and associated evaluations, which is a really challenging issue that is yet to be solved, to the best of our knowledge. In application of variable and model selection, a popular practice is to use the “one standard error rule”. However, its validity hinges crucially on the goodness of the standard error formula used in the approach. Our aim in this work is to study the standard error estimation issue and its consequences on variable selection and regression estimation.

In this section, we provide a background on tuning parameter selection for high-dimensional regression.

1.1. Regularization Methods

Regularization methods are now widely used to tackle the problem of curse of dimensionality in high-dimensional regression analysis. By imposing a penalty on the complexity of the solution, it can solve an ill-posed problem or prevent overfitting. Examples include Lasso ([4]) and Ridge regression ([5,6]), which add an

l_{1}

and

l_{2}

penalty on the coefficient estimates, respectively, to the usual least-squares fit objective function.

Specifically, given a response vector

y \in R^{n}

, a matrix

X \in R^{n \times p}

of predictor variables and a vector of intercept

α

, the Lasso estimate is defined as

{\hat{β}}^{l a s s o} = {argmin}_{β \in R^{p}} {‖ y - α - X β ‖}_{2}^{2} + λ \sum_{j = 1}^{p} ∣ β_{j} ∣,

(1)

and the ridge estimate is defined as

{\hat{β}}^{r i d g e} = {argmin}_{β \in R^{p}} {‖ y - α - X β ‖}_{2}^{2} + λ \sum_{j = 1}^{p} β_{j}^{2} .

(2)

The tuning parameter

λ

controls the strength of the penalty. As is well known, the nature of the

ℓ_{1}

penalty causes some coefficients to be shrunken to zero exactly while the

ℓ_{2}

penalty can shrink all of the coefficients towards zero, but it typically does not set any of them exactly to zero. Thus, an advantage of Lasso is that it makes the results simpler and more interpretable. In this work, ridge regression will be considered when regression estimation (i.e., the focus is on the accurate estimation of the regression function) or prediction is the goal.

1.2. Tuning Parameter Selection

The regularization parameter, or the tuning parameter, determines the extent of penalization. A traditional way to select a model is to use an information criterion such as Akaike information criterion (AIC) ([7,8]), Bayesian information criterion (BIC) and its extension ([9,10,11]). A more commonly used way for tuning parameter selection is the cross-validation ([12,13,14]). For recent work on cross-validation for consistent model selection for high-dimensional regression, see [15].

The idea of K-fold cross validation is to split the data into K roughly equal-sized parts:

F_{1}, \dots, F_{K}

. For the k-th part where

k = 1, 2, \dots, K

, we consider training on

(x_{i}, y_{i}), i \notin F_{k}

and evaluation on

(x_{i}, y_{i}), i \in F_{k}

. For each value of the tuning parameter

λ \in {λ_{1}, \dots, λ_{m}}

, a set of candidate values, we compute the estimate

{\hat{f_{λ}}}^{(- k)}

on the training set, and record the total error on the validation set:

C V_{k} (λ) = \sum_{i \in F_{k}} {(y_{i} - {\hat{f_{λ}}}^{(- k)} (x_{i}))}^{2} .

(3)

For each tuning parameter value

λ

, we compute the average error over all folds:

C V (λ) = \frac{1}{K} \sum_{k = 1}^{K} C V_{k} (λ) .

(4)

We then choose the value of tuning parameter that minimizes the above average error,

{\hat{λ}}_{m i n} = {argmin}_{λ \in {λ_{1}, \dots, λ_{m}}} C V (λ) .

(5)

1.3. Goal for the One Standard Error Rule

The “one standard error rule” (1se rule) for model selection was first proposed for selecting the right sized tree ([16]). It is also suggested to pick the most parsimonious model within one standard error of the minimum cross validation error ([17]). Such a rule acknowledges the fact that the bias-variance trade-off curve is estimated with error, and hence takes a conservative approach: a smaller model is preferred when the prediction errors are more or less indistinguishable. To estimate the standard deviation of

C V (λ)

at each

λ \in {λ_{1}, \dots, λ_{m}}

, one computes the sample standard deviation of validation errors in

C V_{1} (λ), \dots,

C V_{K} (λ)

S D (λ) = \sqrt{V a r (C V_{1} (λ), \dots, C V_{K} (λ))} .

(6)

Then, we estimate the standard deviation of

C V (λ)

, which is declared the standard error of

C V (λ)

, by

S E (λ) = S D (λ) / \sqrt{K} .

(7)

The organization of the rest of the paper is as follows. Section 2 examines validity of the 1se rule from a theoretical perspective. Section 3 describes the objectives and approaches of our numerical investigations. Section 4 presents the experimental results of the 1se rule in three different aspects. Section 5 applies the 1se rule on two real data sets and examine its performances via data guided simulations. Section 6 concludes the paper. Appendix A gives additional numerical results and Appendix B provides the proof of the main theorem in the paper.

2. Theoretical Result on the One-Standard-Error Rule

2.1. When Is the Standard Error Formula Valid?

The 1se rule makes intuitive sense, but its validity is not clear at all. The issue is that the

C V

errors completed at different folds are not independent, which invalidates the standard error formula for the sample mean of i.i.d. observations. Note that the validity of the standard error is equivalent to the validity of the sample standard deviation expressed in Equation (6). Thus, it suffices to study if and to what extent the usual sample variance properly estimates the targeted theoretical variance.

In this section, we investigate the legitimacy of the standard error formula for a general learning procedure. We focus on the regression case, but similar results hold for classification as well.

Let

δ

be a learning procedure that produces the regression estimator

{\hat{f}}_{δ} (x; D^{n})

based on the data

D^{n} = {(x_{i}, y_{i})}_{i = 1}^{n}

. Throughout the section, we assume K is given and for simplicity, n is a multiple of K. Define

C V_{k} (δ) = \sum_{i \in F_{k}} {(y_{i} - {\hat{f}}_{δ} (x_{i}; D_{(- k)}^{n}))}^{2}

, where

D_{(- k)}^{n}

denote the observations in

D^{n}

that do not belong to

F_{k}

and

k = 1, 2, \dots, K

The

S E

of the

C V

errors of the learning procedure

δ

is defined as

S E (δ) = \sqrt{\frac{\frac{1}{K - 1} \sum_{k = 1}^{K} {(C V_{k} (δ) - \frac{1}{K} \sum_{k = 1}^{K} C V_{k} (δ))}^{2}}{K}} .

To study the

S E

formula, given that K is fixed, it is basically equivalent to studying the sample variance formula, i.e.,

S_{n} (δ) = \frac{1}{K - 1} \sum_{k = 1}^{K} {(C V_{k} (δ) - \frac{1}{K} \sum_{k = 1}^{K} C V_{k} (δ))}^{2} .

The validity of these formulas (in an asymptotic sense) hinges on that the correlation between

C V_{k_{1}} (δ)

and

C V_{k_{2}} (δ)

for

k_{1} \neq k_{2}

approaches zero as

n \to \infty

Definition 1.

The sample variance

S_{n} (δ)

is said to be asymptotically relatively unbiased (ARU) if

\frac{E S_{n} (δ) - V a r (C V_{k} (δ))}{V a r (C V_{k} (δ))} \to 0

n \to \infty

Clearly, if the property above holds, the

S E

formula is properly justified in terms of the relative estimation bias being negligible. Then, the 1se rule is sensible.

For

p > 0

, define the

L_{p}

-norm

| | f | | = {(\int f^{p} (x) P_{X} (d x))}^{1 / p},

where

P_{X}

denotes the probability distribution of X from which

X_{1}, X_{2}, \dots, X_{n}

are drawn. Let

{\hat{f}}_{δ, n}

denote the estimator

{\hat{f}}_{δ} (x; D^{n})

Theorem 1.

E (| | f - {\hat{f}}_{δ, n} {| |}_{4}^{4}) = o (\frac{1}{n})

, then

S_{n} (δ)

is asymptotically relatively unbiased.

Remark 1.

The condition

E (| | f - {\hat{f}}_{δ, n} {| |}_{4}^{4}) = o (\frac{1}{n})

holds if the regression procedure is based on a parametric model under regularity conditions. For instance, if we consider a linear model with

p_{0}

terms, then under sensible conditions (e.g., [15], p. 111),

E (| | f - {\hat{f}}_{δ, n} {| |}^{4}) = O (\frac{p_{0}^{2}}{n^{2}})

and consequently as long as

p_{0} = o (\sqrt{n})

, the variance estimator

S_{n} (δ)

is asymptotically relatively unbiased.

Remark 2.

The above result suggests that when relatively simple models are considered, the 1se rule may be reasonably applicable. In the same spirit, when sparsity oriented variable selection methods are considered as candidate learning methods, if the most relevant selected models are quite sparse, the 1se rule may be proper.

2.2. An Illustrative Example

Our theorem shows that the standard error formula is asymptotically valid when the correlation between the

C V

errors at different folds is relatively small, which is related to the rate of convergence of the regression procedure. In this section, we illustrate how the dependence affects the validity of the standard error formula in a concrete example. For ease in highlighting the dependence among the

C V

errors at different folds, we consider a really simple setting, where the true regression function is a constant. The simplicity allows us to calculate the variance and covariance of the

C V

errors explicitly and then see the essence. The insight gained, however, is more broadly applicable.

The true data generating model is

Y = θ + ϵ

, where

θ \in R

is unknown. Given data

D^{n}

, for an integer r as a factor of n, the regression procedure considers only the first, r-th,

(2 r)

-th, ...,

(n / r)

-th observations and take the corresponding sample mean of Y to estimate

f (x) = 0

Suppose

l = \frac{n}{K r}

is an integer. Then, when applying the K-fold

C V

, the training sample size is

(K - 1) \frac{n}{K}

but only

(K - 1) l

observations are actually used to estimate f. Let

W_{i}

W_{j}

denote the

C V

errors on the i-th and j-th fold,

1 \leq i, j \leq K

. Then, it can be shown that

V a r (W_{i}) ≍ n (1 + \frac{n}{l^{2}})

(8)

and

C o v (W_{i}, W_{j}) ≍ \frac{n^{2}}{l^{2}}

(9)

for

i \neq j

, where ≍ denotes the same order.

Therefore, the standard error formula is ARU if and only if

\frac{l}{\sqrt{n}} \to \infty

In this illustrating example, when l is small (recall that only l observations are used in each fold for training), the regression estimator obtained from excluding one fold respectively are highly dependent, which makes the sample variance formula problematic. In contrast, when l is large, e.g.,

l = \frac{n}{K}

(i.e., all the observations are used), the dependence of the

C V

errors at different folds becomes negligible compared to that of the randomness in the response. In the extreme case of

l = \frac{n}{K}

, we have

ρ (W_{i}, W_{j}) ≍ \frac{1}{n}

(10)

for

i \neq j

Note also that the choice of l determines the convergence rate of the regression estimator, i.e.,

E (| | f - {\hat{f}}_{δ, n} | |_{2}^{2})

is of order

1 / l

. Obviously, with large l, the estimator converges fast, which dilutes the dependence of the

C V

errors and helps the sample variance formula to be more reliable.

3. Numerical Research Objectives and Setup

One standard error rule is intended to improve the usual cross-validation method stated in Section 1.2, where

λ

is chosen to minimize the CV error. However, to date, no empirical or theoretical studies were conducted to justify the application of the one standard error rule in favor of parsimony, to the best of our knowledge. There are some interesting open questions that deserve investigations.

(1): Does the 1se itself provide a good estimation of the standard deviation of the cross validation error $C V (λ)$ as intended?
(2): Does the model selected by the 1se rule (model with $λ_{1 s e}$ ) typically outperform the model selected by minimizing the CV error (model with $λ_{m i n}$ ) in variable selection?
(3): What if estimating the regression function or prediction is the goal?

In this paper, we study the estimation accuracy of the one standard error rule and its application in regression estimation and variable selection under various regression scenarios, including both simulations and real data examples.

3.1. Simulation Settings

Our numerical investigations are done in a regression framework where our simulated response is generated by linear models. Both dependent and independent predictors are considered. The data generating process is shown as follows.

Denote the sample size as n and the total number of predictors as p. Firstly, the predictors are generated by normal distributions. In the independent case, the vector of the explanatory variables

(X_{1}, X_{2}, \dots \dots, X_{p})

has i.i.d.

N (0, 1)

components and in the dependent case, the vector follows

N (0, Σ)

distribution, where

Σ

takes either

Σ_{1}

Σ_{2}

as follows

Σ_{1} = (\begin{matrix} 1 & ρ & \dots & ρ^{n - 1} \\ ρ & 1 & \dots & ρ^{n - 2} \\ ⋮ & ⋮ & ⋮ \\ ρ^{n - 1} & ρ^{n - 2} & \dots & 1 \end{matrix})

Σ_{2} = (\begin{matrix} 1 & ρ^{} & \dots & ρ^{} \\ ρ^{} & 1 & \dots & ρ^{} \\ ⋮ & ⋮ & ⋮ \\ ρ^{} & ρ^{} & \dots & 1 \end{matrix}) .

where

0 < ρ < 1

. Then, the response vector

y \in R^{n}

is generated by

y = X β + ϵ,

(11)

where

β

is the chosen coefficient vector (to be specified later) and the random error vector

ϵ

has mean 0 and covariance matrix

σ^{2} I_{n \times n}

for some

σ^{2} > 0

3.2. Regression Estimation

For estimation of the regression function, we intend to minimize the estimation risk. Given the data, we consider a regression method and its loss. For Lasso, we define

{\hat{β}}_{m i n} = {argmin}_{β \in R^{p}} \sum_{i = 1}^{n} {(y_{i} - α_{i} - x_{i}^{'} β)}^{2} + λ \sum_{j = 1}^{p} ∣ β_{j} ∣,

(12)

where

λ

is chosen by the regular CV, i.e.,

λ = λ_{m i n}

. We let

{\hat{β}}_{1 s e}

be the parameter estimate obtained from the simplest model whose CV error is within the one standard error of the minimum CV error, i.e.,

λ = λ_{1 s e}

. Given an estimate

\hat{β}

, the loss for estimating the regression function is

L (\hat{β}, β^{★}) = E {(x^{'} \hat{β} - x^{'} β^{★})}^{2},

(13)

where

β^{★}

is the true parameter vector and the expectation is taken with respect to the randomness of

x

that follows the same distribution as in the data generation. The loss can be approximated by

L (\hat{β}, β^{★}) \approx \frac{1}{J} \sum_{j = 1}^{J} {(z_{j}^{'} \hat{β} - z_{j}^{'} β^{★})}^{2},

(14)

where

z_{j}

j = 1, \dots, J

are i.i.d. with the same distribution as

x

in the data generation and J is large (In the following numerical exercise, we set

J = 500

). In the case of using real data to compare

λ_{1 s e}

and

λ_{m i n}

, clearly, it is not possible to compute the estimation loss as described above. We will use cross-validation and data guided simulations (DGS) instead for comparison of

λ_{1 s e}

and

λ_{m i n}

. See Section 4 for details.

3.3. Variable Selection

Parameter value specification can substantially affect the accuracy of the variable selection by Lasso. In our simulation study, we investigate the performances of the model with

λ_{m i n}

and that with

λ_{1 s e}

by comparing their probabilities of accurate variable selection. When considering real data, since we do not know the true data generation process, we employ data guided simulations to provide more insight. Let

S^{★}

denote the set of variables in the true model and

\hat{S}

be the set of selected variables by a method. In the DGS over real data sets, we define the symmetric difference

S^{★} ▽ \hat{S}

to be the variables that are either in

S^{★}

but not in

\hat{S}

or in

\hat{S}

but not in

S^{★}

. We use the size of the symmetric difference to measure the performance of the variable selection method. Performance of the model with

λ_{m i n}

and the model with

λ_{1 s e}

in terms of variable selection will be compared this way for the data examples, where some plausible models based on the data are used to generate data.

Organization of the numerical results is as follows. Section 4 presents the experimental results of the 1se rule in its estimation accuracy as well as its performances in regression estimation and variable selection, respectively. Section 5 applies the 1se rule over two real data sets and investigates its performance relative to the regular CV via DGS. Some additional figures and tables are in Appendix A. All of R codes used in this paper are available upon request.

4. Simulation Results

4.1. Is 1se on Target of Estimating the Standard Deviation?

This subsection presents a simulation study of the estimation accuracy of the 1se with respect to different parameter settings in three sampling settings, by varying value of the correlation (

ρ

) in the design matrix, error variance

σ^{2}

, sample size n, and the number of cross validation fold K.

4.1.1. Data Generating Processes (DGP)

We consider the following three sampling settings through which we first simulate the dependent (

Σ_{1}

with

ρ = 0.9

) or independent design matrix

X \in R^{n \times p}

, and then generate the response

y = X β + ϵ

with either decaying or constant coefficient vector

β

. Recall that p is the total number of predictors and let q denote the number of nonzero coefficients.

DGP 1: Y is generated by 5 predictors (out of 30, $q = 5$ , $p = 30$ ).
DGP 2: Y is generated by 20 predictors (out of 30, $q = 20$ , $p = 30$ ).
DGP 3: Y is generated by 20 predictors (out of 1000, $q = 20$ , $p = 1000$ ).

For the coefficient $β$ =

(β_{1}, \dots, β_{q}, β_{q + 1}, \dots, β_{p})

, it equals

(1, \frac{1}{2}, \dots, \frac{1}{q}, 0, \dots, 0)

for decaying coefficient case, and it equals

(\underset{q}{\underset{︸}{1, 1, \dots 1}}, \underset{p - q}{\underset{︸}{0, \dots, 0}})

for constant coefficient case. In these three sampling settings,

ϵ

are i.i.d.

N (0, σ^{2})

. Therefore, the first and third configuration yield sparse models and the third configuration also represents the case of high-dimensional data.

4.1.2. Procedure

We apply the Lasso and Ridge procedures on the simulated data sets with tuning parameter selection by the two cross validation versions. The simulation results obtained with various parameter choices are summarized in Table A1, Table A2, Table A3 and Table A4 in the Appendix A. The simulation choices are sample size

n \in {100, 1, 000}

, the error variance

σ^{2} \in {0.01, 1, 4}

and K-fold cross validation:

K \in {2, 5, 10, 20}

(1): Do the j-th simulation ( $j = 1, \dots, N$ ). Simulate a data set of sample size n given the parameter values.
(2): K-fold cross validation. Randomly split the data set into K equal-sized parts: $F_{1}, F_{2}, \dots, F_{K}$ and record the cross validation error for each part $C V_{k} (λ)$ where $k = 1, 2, \dots, K$ .
(3): Calculate the total cross validation error, $C V . e r r_{j} (λ)$ , by taking the mean of these K numbers. Find $λ_{m i n}$ that minimizes the total CV error. Now, we take $λ = λ_{m i n}$ throughout the remaining steps. Calculate the standard error of cross validation error $S E_{j} (λ)$ by one standard error rule:

$S E_{j} (C V . e r r o r) = \sqrt{V a r (C V_{1, j}, C V_{2, j}, \dots, C V_{K, j}) / K} .$
(4): Repeat step 1 to step 3 for N times ( $N = 100$ ). Calculate the standard deviation of the cross validation error $C V . e r r_{j} (λ)$ , and call it $S D (C V . e r r o r)$ , and the mean standard error of cross validation error (as used for 1se rule), $S E (C V . e r r o r) = \frac{1}{N} \sum S E_{j} (C V . e r r o r)$ .
(5): Performance assessment. Calculate the ratio of (the claimed) standard error over (the simulated true) standard deviation:

$\frac{S E (C V . e r r o r)}{S D (C V . e r r o r)} .$

If the 1se works well, the ratio should be close to 1. However, the results in the Table A1, Table A2, Table A3 and Table A4 show that, in some cases, the standard deviation is even more than 100% overestimated or more than 50% underestimated. The histograms of the ratios (Figure A1 and Figure A2) also show that the standard error estimates are not close to the targets. Not surprisingly, 2-fold is particularly untrustworthy. But for the common choice of 10-fold, even for independent predictors with constant coefficient, the standard error estimate by 1se can be 40% over or under the target. In addition, there is no obvious pattern on when the 1se is performing well or not. Overall, this simulation study provides evidence that the 1se fails to provide good estimation as intended (see [18] for more discussions on the standard error issues).

4.2. Is 1se Rule Better for Regression Estimation?

Since the standard error estimation in 1se is not the end goal, the poor performance of 1se in the previous subsection does not necessary mean that the 1se rule is poor for regression estimation and variable selection. In this part, we study the performance of model

λ_{1 s e}

in regression estimation and compare it with that of model

λ_{m i n}

. A large number of samples are drawn so that we know in which cases the 1se rule can consistently outperforms the regular CV in regression estimation.

By simulations, we compare

λ_{1 s e}

and

λ_{m i n}

in term of accuracy in estimating the regression function. We consider several factors: nonzero coefficients

β_{(1)}

, number of zero coefficients

p - q

, sample size n, standard deviation of the noise

σ

and the dependency of the predictors.

4.2.1. Model Specification

We consider the data generation model:

y = X β + ϵ

, where

ϵ \sim N (0, σ^{2})

and

β = (β_{(1)}, β_{(2)}),

with

β_{(1)} = (β_{1}, β_{2}

, \dots, β_{q})

being the nonzero coefficients and

β_{(2)} = (β_{q + 1}, β_{q + 2}

, \dots, β_{p})

being the zero coefficients. The vector of explanatory variables is generated by the multivariate normal distribution,

X \sim M V (0, Σ_{i}),

where

i = 1, 2

4.2.2. Procedure

Simulate a fixed validation set of 500 observations of the predictors for estimation of the loss.
Each time randomly simulate a training set of n observations from the same data generating process.
Apply the 1se rule (10-fold cross validation) over the training set and record the selected model: model with $λ_{1 s e}$ and the model with $λ_{m i n}$ . Calculate the estimation losses of these two models: $\frac{1}{500} \sum_{i = 1}^{500} {(z_{i}^{'} \hat{β} - z_{i}^{'} β^{★})}^{2}$ , where $\hat{β}$ is based on $λ_{1 s e}$ or $λ_{m i n}$ respectively, and $z_{i}$ are independently generated from the same distribution used to generate X.
Repeat this process for M times ( $M = 500$ ) and calculate the fraction that model with $λ_{m i n}$ has a smaller loss.

In our simulations, we make various comparisons: small or large nonzero coefficients; small or large error variance; small or large sample size; low- or high-dimensional data; independent or dependent design matrix. The parameter values used in our study are displayed in Table 1 below. The results are reported in Table A5 and Table A6. Note that when

λ_{m i n}

and

λ_{1 s e}

give identical selection over the

M = 500

runs, the entry in the table is NaN.

The simulation results show that, in general, for the design matrix with either AR(1) correlation or constant correlation, the model with

λ_{m i n}

tends to provide better regression estimation (with smaller estimation errors), especially when the data is generated with relatively large coefficient. In addition, we consider two more settings below.

4.2.3. Case 1: Large $p / q$

In this case, we consider high-dimensional data with AR(1) covariance design matrix,

Σ_{1}

, where

q = 10

and

p = 1000

. We set

σ^{2} = 1

n \in {100, 1000}

ρ \in {0.01, 0.9}

, constant nonzero coefficients with value

β \in {1, 6, 11, 16, 21}

. We report the probability that model

λ_{m i n}

estimates better. Clearly, in this case, the model

λ_{m i n}

significantly outperforms the model

λ_{1 s e}

for a range of

β

values at both sample size (See Figure 1).

4.2.4. Case 2: Small $p / q$

In this case, we consider low-dimensional data with AR(1) covariance design matrix,

Σ_{1}

, where

q = 10

and

p - q \in {2, 5, 10}

. We set

σ^{2} = 0.01

n \in {100, 1000}

ρ = 0.1

, the common nonzero coefficient value

β \in {0.2, 0.4, 0.6, 0.8, 1,

1.2}

. We report the probability that model

λ_{m i n}

estimates the regression function better. In this case, model

λ_{m i n}

significantly outperforms model

λ_{1 s e}

as the common coefficient value increases over 1, getting close to or even over 0.8 (See Figure 2).

4.3. Is 1se Rule Better for Variable Selection?

There are two main issues of the application of the 1se rule in variable selection. One is that the accuracy of variable selection by Lasso itself is very sensitive to the parameter values in the data generating process. In various situations (especially when the covariance are complicatedly correlated such as in gene expression type of data), Lasso may perform very poorly in accurate variable selection. The second issue is that

λ_{m i n}

and

λ_{1 s e}

returned by the glmnet function are the same in many cases and then there is no difference in the result of variable selection. The cases where

λ_{m i n} = λ_{1 s e}

are eliminated in the following general study. To examine the real difference between

λ_{m i n}

and

λ_{1 s e}

, we consider the difference of their probabilities of selecting the true model given that they have distinct select outcome.

4.3.1. Procedure

Randomly simulate a data set of sample size n given the parameter values.
Perform variable selection with $λ_{m i n}$ and $λ_{1 s e}$ over the simulated data set, respectively. If the two returned models are the same, then discard the result and go back to step 1. Otherwise, check their variable selection results.
Repeat the above process for M times ( $M = 500$ ) and calculate the fraction of correct variable selection, denoted by $P_{c} (λ_{1 s e})$ and $P_{c} (λ_{m i n})$ , respectively. We report the proportion difference $P_{c} (λ_{1 s e}) - P_{c} (λ_{m i n})$ . Positive means 1se rule is better.

We adopt the same framework and parameter settings (Table 1) as that used for regression estimation. The results show that

λ_{1 s e}

tends to outperform

λ_{m i n}

in most of the cases when the data are generated with AR(1) correlation design matrix (see Table A7 and Table A8). But with constant correlation design matrix, these two methods have no practically difference in variable selection.

Below we consider several special cases where model

λ_{m i n}

still needs to be considered, and although it can only do slightly better when the error variance is large, nonzero coefficients are small and no extra variables exist to interrupt the selection process.

4.3.2. Case 1: Constant Coefficients

In this case, let

q = 5

p - q \in {0, 1, 2, 3, 4}

σ^{2} = 4

n \in {100, 1000}

ρ = 0.01

and the common nonzero

β \in [0.5, 1.3]

with

Σ_{1}

. The result is in Figure 3. In the case of small sample size, say

n = 100

, model

λ_{m i n}

has higher probability of correct variable selection when the constant coefficient lies between (0.5,1). The largest difference appears when

p - q = 0

, not surprisingly. As the common coefficient value becomes larger than 1, the probability increases from negative to positive, which indicates the reverse to be true: model

λ_{1 s e}

now works better. However, in the case of large sample size, say

n = 1000

, model

λ_{1 s e}

significantly outperforms model

λ_{m i n}

unless when

p - q = 0

4.3.3. Case 2: Decaying Coefficients

In this case, let

q = 5

σ^{2} = 4

n = 1000

ρ = 0.5

p - q \in 0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50

, and

β_{i} = \frac{1}{i}

i = 1, 2, \dots, q

with

Σ_{1}

. Figure 4 indicates that, with decaying coefficients lie between (0,1) and large sample size

n = 1000

, model

λ_{m i n}

has higher probability of correct variable selection when the number of zero-coefficient variables

p - q

is small. The largest difference appears when

p - q = 0

, i.e., there is no extra variables provided to interrupt the variable selection process.

4.3.4. Case 3: Hybrid Case

Let

n = 100

σ^{2} = 4

and

ρ \in {0, 0.5}

with

Σ_{1}

. In this data generating process, we set the first ten coefficients to be positive, i.e.,

q = 10

, among which the first nine nonzero coefficients are constant 1, while the last one

β_{10}

increases from 0 to 10, i.e.,

β_{10} \in {0, 0.05, 0.1, 0.15, 0.2, 0.4, 0.6, 0.8, 1, 3, 5, 10}

. The number of zero coefficients

p - q \in {0, 1, 2, 3,

5, 20}

. From Figure 5, when

β_{10}

is very small, i.e.,

β_{10} \in {0.05, 0.2}

and

p - q = 0

λ_{m i n}

performs better. For other situations,

λ_{1 s e}

is significantly better.

5. Data Examples

In this section, the performances of the 1se rule and the regular CV in regression estimation and variable selection are examined over real data sets: Boston Housing Price (an example of

n > p

) and Bardet–Biedl data (an example of

p > n

Boston Housing Price. The data were collected by [19] for the purpose of discovering whether or not clean air influenced the value of houses in Boston. It consists of 506 observations and 14 non constant independent variables. Of these, medv is the response variable while the other 13 variables are possible predictors.

Bardet–Biedl. The gene expression data are from the microarray experiments of mammalianeye tissue samples of [20]. The data consists of 120 rats with 200 gene probes and the expression level of TRIM32 gene. It is interesting to know which gene probes are more able to explain the expression level of TRIM32 gene by regression analysis.

5.1. Regression Estimation

To test which method does better in regression estimation, we perform both cross validation and data guided simulation (DGS) over these two real data sets as follows.

5.1.1. Procedure for Cross Validation

We randomly select $n_{1}$ observations from the data sets as training set and the rest as validation set, $n_{1} \in {100, 200,$ $400}$ for Boston Housing Price and $n_{1} \in {40, 80}$ for Bardet–Biedl data.
Apply K-fold cross validation $K \in {5, 10, 20}$ , for Boston Housing Price and $K \in {5, 10}$ for Bardet–Biedl data over the training set and compute the mean square prediction error of model $λ_{1 s e}$ and model $λ_{m i n}$ over the validation set.
Repeat the above process for 500 times and compute the proportion of better estimation for each method.

The results (Table 2 and Table 3) show that, for all the training size

n_{1}

considered, the proportion of model

λ_{m i n}

doing better is about 52∼57% for Boston Housing Price; 56∼59% for Bardet–Biedl data. The consistent slight advantage of

λ_{m i n}

agrees with our earlier simulation results in supporting

λ_{m i n}

as the better method for regression estimation.

5.1.2. Procedure for DGS

Obtain model $λ_{m i n}$ and model $λ_{1 s e}$ with K-fold cross validation, $K \in {5, 10, 20}$ , over the data and also the coefficient estimates ${\hat{β}}_{m i n}$ , ${\hat{β}}_{1 s e}$ , the standard deviation estimates ${\hat{σ}}_{m i n}$ , ${\hat{σ}}_{1 s e}$ and the estimated responses ${\hat{y}}_{m i n, 1}$ , ${\hat{y}}_{1 s e, 1}$ (fitted values).
Simulate the new responses in two scenarios:
Scenario 1: $y = {\hat{y}}_{m i n, 1} + ϵ$ with $ϵ \sim N (0, {\hat{σ}}_{m i n}^{2} I)$ ,
Scenario 2: $y = {\hat{y}}_{1 s e, 1} + ϵ$ with $ϵ \sim N (0, {\hat{σ}}_{1 s e}^{2} I)$ .
Here, apply Lasso with $λ_{m i n}$ and $λ_{1 s e}$ (using the same K-fold CV) on the new data set (i.e., new response and the original design matrix) and get the new estimated responses ${\hat{y}}_{m i n, 2}$ and ${\hat{y}}_{1 s e, 2}$ for each of the two scenarios.
Calculate the mean square estimation error $\frac{1}{n} \sum_{i = 1}^{n} {({\hat{y}}_{m i n, 2} - {\hat{y}}_{m i n, 1})}^{2}$ and $\frac{1}{n} \sum_{i = 1}^{n} {({\hat{y}}_{1 s e, 2} - {\hat{y}}_{m i n, 1})}^{2}$ for scenario 1;
$\frac{1}{n} \sum_{i = 1}^{n} {({\hat{y}}_{m i n, 2} - {\hat{y}}_{1 s e, 1})}^{2}$ and $\frac{1}{n} \sum_{i = 1}^{n} {({\hat{y}}_{1 s e, 2} - {\hat{y}}_{1 s e, 1})}^{2}$ for scenario 2.
Repeat the above resampling process for 500 times and compute the proportion of better estimation for each method.

The estimation results (Table 4) show that, overall, the proportion of better estimation for the model

λ_{m i n}

is between 46.8∼52.8% and for the model

λ_{1 s e}

, it is between 47.2∼53.2%. Therefore, based on the DGS, in term of regression estimation accuracy, there is not much difference before

λ_{m i n}

and

λ_{1 s e}

: the observed properties are not significantly different from 50% at 0.05 level.

5.2. Variable Selection: DGS

To test which model does better in variable selection, we perform the DGS over the above real data sets. Since the exact variable selection might not be reached, especially when applied to the high-dimensional data, we use symmetric difference

∣ S^{★} ▽ \hat{S} ∣

to measure their performance, which is the number of variables either in

S^{★}

but not in

\hat{S}

or in

\hat{S}

but not in

S^{★}

. Below is the algorithm for the DGS.

Procedure

Apply the 10-fold (default) cross validation over real data set and select the true set of variables $S^{★}$ : $S_{m i n}^{★}$ or $S_{1 s e}^{★}$ , either by the model with $λ_{m i n}$ or by the model with $λ_{1 s e}$ .
Do least square estimation by regressing the response on the true set of variables and obtain the estimated response $\hat{y}$ and the residual standard error $\hat{σ}$ .
Simulate the new response $y_{n e w}$ by adding the error term randomly generated from $N (0, {\hat{σ}}^{2})$ to $\hat{y}$ . Apply K-fold cross validation, $K \in {5, 10, 20}$ , over the simulated data set (i.e., $y_{n e w}$ and the original design matrix) and select set of variables $\hat{S}$ : ${\hat{S}}_{m i n}$ or ${\hat{S}}_{1 s e}$ , either by model $λ_{m i n}$ or by model $λ_{1 s e}$ . Repeat this process for 500 times.
Calculate the symmetric difference $∣ S^{★} ▽ \hat{S} ∣$ .

The distributions of the symmetric difference are shown in the appendix. (Figure A3 and Figure A4 for Boston Housing Price; Figure A5 and Figure A6 for Bardet–Biedl data) The mean of

∣ \hat{S} ∣

is reported in Table 5 and Table 6. For Boston Housing Price, the size of the candidate variables is 13 while model

λ_{m i n}

selects 11 variables and model

λ_{1 s e}

selects 8 variables on average, i.e., the mean of

∣ S_{m i n}^{★} ∣

is 11 and the mean of

∣ S_{1 s e}^{★} ∣

is 8. For Bardet–Biedl data, the size of the candidate variables is 200, while model

λ_{m i n}

selects 21 variables and model

λ_{1 s e}

selects 19 variables on average, i.e., the mean of

∣ S_{m i n}^{★} ∣

is 21 and the mean of

∣ S_{1 s e}^{★} ∣

is 19.

On average, the 1se rule does better in variable selection since the mean of

∣ {\hat{S}}_{1 s e} ∣

is more closed to

∣ S^{★} ∣

. In the case of

n > p

, if we want to have a perfect or near perfect variable selection result (i.e.,

∣ S^{★} ▽ \hat{S} ∣ \leq 1

), then the 1se rule is a good choice, despite the fact that its selection result is more unstable: the range of

∣ S^{★} ▽ {\hat{S}}_{1 s e} ∣

is larger than that of

∣ S^{★} ▽ {\hat{S}}_{m i n} ∣

. In the case of

n < p

, both models fail to provide perfect selection result. But overall the 1se rule is more accurate and stable given that

∣ S^{★} ▽ {\hat{S}}_{1 s e} ∣

has smaller mean and range.

6. Conclusions

The one standard error rule is proposed to pick the most parsimonious model within one standard error of the minimum. Despite the fact that it was widely used as a conservative and alternative approach for cross validation, there is no evidence confirming that the 1se rule can consistently outperform the regular CV.

Our theoretical result shows that the standard error formula is asymptotically valid when the regression procedure converges relatively fast to the true regression function. Our illustrative example also shows that when the regression estimator converges slowly, the CV errors on the different evaluation folds may be highly dependent, which unfortunately makes the usual sample variance (or sample standard deviation) formula problematic. In such a case, the use of the 1se rule may easily lead to a choice of inferior candidates.

Our numerical results offer finite-sample understandings. First of all, the simulation study also casts doubt on the estimation accuracy of the 1se rule in terms of estimating the standard error of the minimum cross validation error. In some cases, the 1se is likely to have 100% overestimation or underestimation, depending on the number of cross validation fold and the dependency among the predictors. More specifically, for example, the 1se tends to have overestimation when used with 2-fold cross validation in case of dependent predictors; but it tends to have underestimation when used with 20-fold cross validation in case of independent predictors.

Secondly, the performances of the 1se rule and the regular CV in regression estimation and variable selection are compared via both simulated and real data sets. In general, on the one hand, the 1se rule often performs better in variable selection for sparse modeling; on the other hand, it often does worse in regression estimation.

While our theoretical and numerical results clearly challenge indiscriminate uses of the 1se rule in general model selection, its parsimony nature may still be appealing when a sparse model is highly desirable. For the real data sets, for regression estimation, the regular CV does better than the 1se rule but the difference is small. There is almost no difference in the estimation accuracy between these two methods when regression estimation is assessed by the DGS. Overall, considering the model simplicity and interpretability, the 1se rule would be a better choice here.

Future studies are encouraged to provide more theoretical understandings on the 1se rule. Also, it may be studied numerically in more complex settings (e.g., when used to compare several nonlinear/nonparametric classification methods) than those in this work. In a larger context, since it is well known that the penalty methods for model selection are often unstable (e.g., [1,21,22]), alternative methods to strengthen reliability and reproductivity of variable selection may be applied (see reference in the above mentioned papers and [23]).

Author Contributions

Conceptualization, Y.Y.; methodology, Y.C. and Y.Y.; numerical results, Y.C.; writing—original draft preparation, Y.C.; writing—review and editing, Y.Y. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Additional Numerical Results

Figure A1. Histogram of Table A1 and Table A2.

Figure A2. Histogram of Table A3 and Table A4.

Table A1. Mean ratio of claimed SE and true SD: dependent predictors (

Σ_{1}

) and decaying coefficients.

Table A1. Mean ratio of claimed SE and true SD: dependent predictors (

Σ_{1}

) and decaying coefficients.

			DGP1				DGP2				DGP3
			2-fold	5-fold	10-fold	20-fold	2-fold	5-fold	10-fold	20-fold	2-fold	5-fold	10-fold	20-fold
$n = 100$	$σ^{2} = 0.01$	Lasso	1.489	1.305	1.275	1.271	1.068	1.185	1.183	1.127	1.015	1.217	1.266	1.222
	$σ^{2} = 0.01$	Ridge	1.187	1.104	0.975	0.924	1.370	1.171	1.082	1.031	0.918	1.190	1.252	1.194
	$σ^{2} = 1$	Lasso	1.759	1.376	1.145	1.088	1.247	1.181	1.148	1.134	0.886	1.094	0.906	0.728
	$σ^{2} = 1$	Ridge	1.499	1.350	1.231	1.129	1.512	1.310	1.202	1.102	1.003	1.092	1.093	1.049
	$σ^{2} = 4$	Lasso	1.260	1.162	1.067	1.071	1.333	1.202	1.164	1.119	1.121	0.580	0.319	0.229
	$σ^{2} = 4$	Ridge	1.540	1.332	1.217	1.108	1.553	1.360	1.252	1.142	1.080	0.969	0.968	0.949
$n = 1000$	$σ^{2} = 0.01$	Lasso	1.221	1.226	1.123	1.134	1.876	1.100	1.028	0.997	1.519	1.166	1.101	1.111
	$σ^{2} = 0.01$	Ridge	1.725	1.011	0.880	0.853	1.490	0.882	0.780	0.777	1.385	1.299	1.254	1.251
	$σ^{2} = 1$	Lasso	1.960	1.233	1.068	1.031	1.725	1.219	1.073	0.959	1.853	0.889	0.711	0.637
	$σ^{2} = 1$	Ridge	2.173	1.243	1.060	1.029	2.069	1.188	1.010	0.978	1.346	1.224	1.183	1.177
	$σ^{2} = 4$	Lasso	1.650	1.031	0.955	0.873	2.122	1.184	1.003	0.880	1.868	0.831	0.668	0.598
	$σ^{2} = 4$	Ridge	2.139	1.218	1.032	0.996	2.154	1.244	1.057	1.020	1.180	1.104	1.085	1.079

Table A2. Mean ratio of claimed SE and true SD: independent predictors and decaying coefficients.

			DGP1				DGP2				DGP3
			2-fold	5-fold	10-fold	20-fold	2-fold	5-fold	10-fold	20-fold	2-fold	5-fold	10-fold	20-fold
$n = 100$	$σ^{2} = 0.01$	Lasso	1.275	1.225	1.130	1.077	1.029	1.087	1.110	1.128	1.221	0.952	0.581	0.379
	$σ^{2} = 0.01$	Ridge	0.973	0.987	0.918	0.909	1.035	1.012	0.964	0.905	1.152	1.130	1.143	1.108
	$σ^{2} = 1$	Lasso	1.532	1.197	1.055	0.999	1.121	1.054	1.132	1.030	1.243	0.978	0.651	0.459
	$σ^{2} = 1$	Ridge	1.039	0.989	0.936	0.924	1.083	1.043	0.977	0.946	1.049	0.986	1.004	0.967
	$σ^{2} = 4$	Lasso	1.569	1.241	1.074	1.026	1.109	1.125	1.061	1.041	1.236	1.033	0.767	0.584
	$σ^{2} = 4$	Ridge	1.027	0.990	0.938	0.933	1.050	1.030	0.963	0.946	1.006	0.951	0.963	0.930
$n = 1000$	$σ^{2} = 0.01$	Lasso	1.126	0.914	0.748	0.668	1.499	1.046	0.864	0.820	1.855	0.961	0.789	0.741
	$σ^{2} = 0.01$	Ridge	1.509	1.108	0.964	0.909	1.677	1.007	0.900	0.854	1.203	1.127	1.100	1.091
	$σ^{2} = 1$	Lasso	1.100	0.749	0.641	0.597	1.624	1.181	0.984	0.888	2.065	1.026	0.823	0.747
	$σ^{2} = 1$	Ridge	1.264	1.037	0.975	0.959	1.327	0.994	0.944	0.926	1.003	0.962	0.962	0.964
	$σ^{2} = 4$	Lasso	1.864	1.193	0.964	0.879	1.866	0.965	0.797	0.710	2.102	1.042	0.836	0.756
	$σ^{2} = 4$	Ridge	1.145	1.004	0.963	0.960	1.131	0.975	0.944	0.938	0.955	0.935	0.937	0.941

Table A3. Mean ratio of claimed SE and true SD: dependent predictors (

Σ_{1}

) and constant coefficients.

Table A3. Mean ratio of claimed SE and true SD: dependent predictors (

Σ_{1}

) and constant coefficients.

			DGP1				DGP2				DGP3
			2-fold	5-fold	10-fold	20-fold	2-fold	5-fold	10-fold	20-fold	2-fold	5-fold	10-fold	20-fold
$n = 100$	$σ^{2} = 0.01$	Lasso	1.452	1.346	1.322	1.307	0.958	1.132	1.140	1.103	0.953	1.155	1.214	1.161
	$σ^{2} = 0.01$	Ridge	1.211	1.279	1.068	1.047	1.573	1.391	1.246	1.184	1.046	1.370	1.399	1.356
	$σ^{2} = 1$	Lasso	1.595	1.284	1.170	1.134	0.955	1.186	1.195	1.136	0.947	1.178	1.231	1.203
	$σ^{2} = 1$	Ridge	1.307	1.280	1.144	1.074	1.783	1.428	1.384	1.290	1.068	1.397	1.440	1.380
	$σ^{2} = 4$	Lasso	1.529	1.415	1.176	1.070	0.998	1.280	1.260	1.254	1.039	1.312	1.290	1.270
	$σ^{2} = 4$	Ridge	1.432	1.340	1.221	1.130	1.566	1.312	1.208	1.111	1.092	1.422	1.471	1.389
$n = 1000$	$σ^{2} = 0.01$	Lasso	1.193	1.187	1.110	1.119	1.084	0.804	0.773	0.744	1.492	1.128	1.122	1.079
	$σ^{2} = 0.01$	Ridge	2.087	1.123	0.970	0.912	1.742	1.448	1.423	1.301	1.726	1.596	1.509	1.482
	$σ^{2} = 1$	Lasso	1.995	1.469	1.330	1.292	1.389	1.035	0.950	0.948	1.573	1.225	1.154	1.159
	$σ^{2} = 1$	Ridge	2.017	1.196	1.029	1.011	2.315	1.448	1.200	1.088	1.744	1.599	1.505	1.486
	$σ^{2} = 4$	Lasso	2.200	1.023	0.883	0.806	2.128	1.375	1.162	1.127	2.144	1.385	1.235	1.204
	$σ^{2} = 4$	Ridge	2.163	1.260	1.078	1.054	1.991	1.261	1.050	1.008	1.783	1.612	1.500	1.490

Table A4. Mean ratio of claimed SE and true SD: independent predictors and constant coefficients.

			DGP1				DGP2				DGP3
			2-fold	5-fold	10-fold	20-fold	2-fold	5-fold	10-fold	20-fold	2-fold	5-fold	10-fold	20-fold
$n = 100$	$σ^{2} = 0.01$	Lasso	1.180	1.248	1.178	1.158	0.921	1.082	1.160	1.109	1.197	1.172	0.599	0.386
	$σ^{2} = 0.01$	Ridge	0.918	1.149	1.099	1.104	1.121	1.089	1.009	0.851	1.103	1.112	1.096	1.050
	$σ^{2} = 1$	Lasso	1.311	1.241	1.133	1.075	0.930	1.110	1.195	1.129	1.371	1.166	0.594	0.380
	$σ^{2} = 1$	Ridge	0.953	0.980	0.898	0.886	1.116	0.997	0.858	0.740	1.122	1.116	1.096	1.054
	$σ^{2} = 4$	Lasso	1.366	1.202	1.087	1.057	0.966	1.104	1.151	1.123	1.387	1.155	0.597	0.382
	$σ^{2} = 4$	Ridge	1.018	0.985	0.917	0.906	1.135	1.016	0.896	0.799	1.159	1.107	1.092	1.052
$n = 1000$	$σ^{2} = 0.01$	Lasso	0.939	0.978	0.842	0.763	1.149	1.030	1.093	1.003	1.068	1.050	1.057	1.067
	$σ^{2} = 0.01$	Ridge	1.475	0.969	0.768	0.719	1.101	1.045	0.942	0.896	1.543	1.446	1.412	1.396
	$σ^{2} = 1$	Lasso	1.311	1.003	0.833	0.788	0.991	1.001	1.065	0.954	1.181	1.027	1.063	1.017
	$σ^{2} = 1$	Ridge	1.511	1.086	0.924	0.874	1.362	0.572	0.460	0.433	1.370	1.185	1.154	1.142
	$σ^{2} = 4$	Lasso	0.996	0.824	0.650	0.631	1.317	0.971	0.909	0.798	1.293	1.000	1.014	0.957
	$σ^{2} = 4$	Ridge	1.405	1.099	0.974	0.934	1.552	0.793	0.679	0.652	1.303	1.157	1.136	1.131

Table A5. Regression estimation: probability of

λ_{m i n}

doing better, AR(1) correlation.

Table A5. Regression estimation: probability of

λ_{m i n}

doing better, AR(1) correlation.

			Low-Dimension				High-Dimension
			$β = 0.1$	$β = 1$	$β = 3$	Decay $β$	$β = 0.1$	$β = 1$	$β = 3$	Decay $β$
$n = 100$	$ρ = 0$	$σ^{2} = 0.01$	0.52	0.62	1.00	0.52	0.46	0.64	0.84	0.54
	$ρ = 0$	$σ^{2} = 4$	0.42	0.52	0.52	0.54	0.58	0.50	0.50	0.52
	$ρ = 0.5$	$σ^{2} = 0.01$	0.54	0.83	0.92	0.52	0.58	0.70	0.90	0.58
	$ρ = 0.5$	$σ^{2} = 4$	0.44	0.56	0.54	0.58	0.38	0.60	0.60	0.60
	$ρ = 0.9$	$σ^{2} = 0.01$	0.56	0.88	1.00	0.60	0.54	0.84	0.96	0.54
	$ρ = 0.9$	$σ^{2} = 4$	0.54	0.56	0.56	0.54	0.57	0.52	0.52	0.46
$n = 1000$	$ρ = 0$	$σ^{2} = 0.01$	0.58	0.57	NaN	0.58	0.64	0.71	NaN	0.74
	$ρ = 0$	$σ^{2} = 4$	0.45	0.58	0.56	0.52	0.57	0.66	0.66	0.72
	$ρ = 0.5$	$σ^{2} = 0.01$	0.50	NaN	NaN	0.54	0.68	NaN	NaN	0.68
	$ρ = 0.5$	$σ^{2} = 4$	0.52	0.50	0.52	0.54	0.64	0.68	0.68	0.70
	$ρ = 0.9$	$σ^{2} = 0.01$	0.56	NaN	NaN	0.60	0.66	NaN	NaN	0.70
	$ρ = 0.9$	$σ^{2} = 4$	0.56	0.54	0.58	0.54	0.68	0.66	0.64	0.66

Table A6. Regression estimation: probability of

λ_{m i n}

doing better, constant correlation.

Table A6. Regression estimation: probability of

λ_{m i n}

doing better, constant correlation.

			Low-Dimension				High-Dimension
			$β = 0.1$	$β = 1$	$β = 3$	Decay $β$	$β = 0.1$	$β = 1$	$β = 3$	Decay $β$
$n = 100$	$ρ = 0.5$	$σ^{2} = 0.01$	0.56	0.88	0.89	0.53	0.54	0.84	1.00	0.64
	$ρ = 0.5$	$σ^{2} = 4$	0.53	0.56	0.54	0.52	0.61	0.52	0.56	0.60
	$ρ = 0.9$	$σ^{2} = 0.01$	0.60	1.00	1.00	0.70	0.56	0.82	0.90	0.70
	$ρ = 0.9$	$σ^{2} = 4$	0.58	0.54	0.62	0.62	0.62	0.58	0.64	0.64
$n = 1000$	$ρ = 0.5$	$σ^{2} = 0.01$	0.50	NaN	NaN	0.50	0.74	NaN	NaN	0.78
	$ρ = 0.5$	$σ^{2} = 4$	0.50	0.42	0.50	0.44	0.72	0.74	0.76	0.78
	$ρ = 0.9$	$σ^{2} = 0.01$	0.48	NaN	NaN	0.54	0.74	NaN	NaN	0.78
	$ρ = 0.9$	$σ^{2} = 4$	0.52	0.50	0.52	0.46	0.72	0.74	0.76	0.76

Table A7. Variable selection:

P_{c} (λ_{1 s e}) - P_{c} (λ_{m i n})

, AR(1) correlation.

Table A7. Variable selection:

P_{c} (λ_{1 s e}) - P_{c} (λ_{m i n})

, AR(1) correlation.

			Low-Dimension				High-Dimension
			$β = 0.1$	$β = 1$	$β = 3$	Decay $β$	$β = 0.1$	$β = 1$	$β = 3$	Decay $β$
$n = 100$	$ρ = 0$	$σ^{2} = 0.01$	0.13	0.13	0	0.06	0	0	0.02	NaN
	$ρ = 0$	$σ^{2} = 4$	0	0.12	0.11	0.01	0	0	0	0.90
	$ρ = 0.5$	$σ^{2} = 0.01$	0.45	0	0	0	0.25	0.01	0	NaN
	$ρ = 0.5$	$σ^{2} = 4$	0	0.40	0.43	0.05	0	0.16	0.27	0.87
	$ρ = 0.9$	$σ^{2} = 0.01$	0.40	0	0	0	0.61	0	0	NaN
	$ρ = 0.9$	$σ^{2} = 4$	0	0.14	0.39	0.01	0	0.25	0.62	0.7
$n = 1000$	$ρ = 0$	$σ^{2} = 0.01$	0.72	0.10	NaN	0.70	0.27	0.16	NaN	NaN
	$ρ = 0$	$σ^{2} = 4$	0	0.71	0.70	0	0	0.25	0.25	0.52
	$ρ = 0.5$	$σ^{2} = 0.01$	0.86	NaN	NaN	0.03	0.86	NaN	NaN	NaN
	$ρ = 0.5$	$σ^{2} = 4$	0	0.85	0.85	0	0	0.85	0.85	0.95
	$ρ = 0.9$	$σ^{2} = 0.01$	0.56	NaN	NaN	0	0.76	NaN	NaN	NaN
	$ρ = 0.9$	$σ^{2} = 4$	0	0.67	0.36	0	0	0.77	0.17	0.65

Table A8. Variable selection:

P_{c} (λ_{1 s e}) - P_{c} (λ_{m i n})

: constant correlation.

Table A8. Variable selection:

P_{c} (λ_{1 s e}) - P_{c} (λ_{m i n})

: constant correlation.

			Low-Dimension				High-Dimension
			$β = 0.1$	$β = 1$	$β = 3$	Decay $β$	$β = 0.1$	$β = 1$	$β = 3$	Decay $β$
$n = 100$	$ρ = 0.5$	$σ^{2} = 0.01$	0	0	0.1	0	0	0	0	0
	$ρ = 0.5$	$σ^{2} = 4$	0	0	0	0	0	0	0	0
	$ρ = 0.9$	$σ^{2} = 0.01$	0	0	0	0	0	0	0	0
	$ρ = 0.9$	$σ^{2} = 4$	0	0	0	0	0	0	0	0
$n = 1000$	$ρ = 0.5$	$σ^{2} = 0.01$	0	NaN	NaN	0.01	0	NaN	NaN	0
	$ρ = 0.5$	$σ^{2} = 4$	0	0	0	0	0	0	0	0
	$ρ = 0.9$	$σ^{2} = 0.01$	0	NaN	NaN	0.01	0	NaN	NaN	0
	$ρ = 0.9$	$σ^{2} = 4$	0	0	0	0	0	0	0	0

Figure A3. Boston Housing Price:

∣ S ▽ \hat{S} ∣; S = v a r_{m i n},

(a)

\hat{S} = v a r_{m i n};

(b)

\hat{S} = v a r_{1 s e}

Figure A3. Boston Housing Price:

∣ S ▽ \hat{S} ∣; S = v a r_{m i n},

(a)

\hat{S} = v a r_{m i n};

(b)

\hat{S} = v a r_{1 s e}

Figure A4. Boston Housing Price:

∣ S ▽ \hat{S} ∣; S = v a r_{1 s e},

(a)

\hat{S} = v a r_{m i n};

(b)

\hat{S} = v a r_{1 s e}

Figure A4. Boston Housing Price:

∣ S ▽ \hat{S} ∣; S = v a r_{1 s e},

(a)

\hat{S} = v a r_{m i n};

(b)

\hat{S} = v a r_{1 s e}

Figure A5. Bardet–Biedl data:

∣ S ▽ \hat{S} ∣; S = v a r_{m i n},

(a)

\hat{S} = v a r_{m i n};

(b)

\hat{S} = v a r_{1 s e}

Figure A5. Bardet–Biedl data:

∣ S ▽ \hat{S} ∣; S = v a r_{m i n},

(a)

\hat{S} = v a r_{m i n};

(b)

\hat{S} = v a r_{1 s e}

Figure A6. Bardet–Biedl data:

∣ S ▽ \hat{S} ∣; S = v a r_{1 s e},

(a)

\hat{S} = v a r_{m i n};

(b)

\hat{S} = v a r_{1 s e}

Figure A6. Bardet–Biedl data:

∣ S ▽ \hat{S} ∣; S = v a r_{1 s e},

(a)

\hat{S} = v a r_{m i n};

(b)

\hat{S} = v a r_{1 s e}

Appendix B. Proof of the Main Result

Appendix B.1. A Lemma

Lemma A1.

Let

W_{1}, W_{2}, \dots, W_{K}

have identical variance and common pairwise correlation ρ. Then

E \frac{1}{K - 1} \sum_{i = 1}^{K} {(W_{i} - \bar{W})}^{2} = (1 - ρ) V a r (W_{i}) .

(A1)

Remark A1.

Depending on

ρ = 0

ρ > 0

ρ < 0

, the sample variance is unbiased, under-estimating, or over estimating the variance of

W_{i}

. For a general learning procedure, the correlation between the test errors from different folds can be both positive and negative, which makes the estimated standard error based on the usual sample variance formula under-estimating or over-estimating the true standard error. This is seen in our numerical results.

Remark A2.

Note that for the sample mean

\bar{W}

, we have

V a r (\bar{W}) = V a r (W_{1}) / K + V a r (W_{1}) ρ (K - 1) / K

. Thus for

V a r (W_{1}) / K

to be asymptotically valid (in approximation) for

V a r (\bar{W})

in the sense that their ratio approaches 1, we must have

ρ \to 0

, which is also the key requirement for the sample variance formula to be asymptotically valid (ARU).

Appendix B.2. Proof of Lemma A1

Since the sample variance stays the same with a common shift of the observations, without loss of generality, we assume

E W_{i} = 0

. Then

\begin{matrix} \begin{matrix} S^{2} & = \frac{1}{K - 1} \sum_{i = 1}^{K} {(W_{i} - \bar{W})}^{2} \\ = \frac{1}{K - 1} (\sum_{i = 1}^{K} W_{i}^{2} - K {\bar{W}}^{2}) \\ = \frac{1}{K - 1} (\sum_{i = 1}^{K} W_{i}^{2} - \frac{1}{K} \sum W_{i}^{2} - \frac{2}{K} \sum_{1 \leq i < j \leq K} W_{i} W_{j}) \\ = \frac{1}{K} \sum_{i = 1}^{K} W_{i}^{2} - \frac{2}{K (K - 1)} \sum_{i < j} W_{i} W_{j} . \end{matrix} \end{matrix}

(A2)

Thus we conclude

E S^{2} = (1 - ρ) V a r (W_{i})

Appendix B.3. Proof of Theorem 1

From Lemma A1, we only need to show that the correlation between the test errors on any two folds approaches 0 as

n \to \infty

. In the following, we examine the covariance between two test errors and show it is negligible compared to the variance of the test error under the condition in the theorem.

Let

W_{1}, \dots, W_{K}

denote the

C V

test errors for the k-th fold,

k = 1, 2, \dots, K

. Note

W_{k} = \sum_{j \in F_{k}} {(Y_{j} - {\hat{f}}_{(- F_{k})} (X_{j}))}^{2}

, where

{\hat{f}}_{(- F_{k})}

denote the estimator of f by the learning procedure

δ

based on

D_{(- F_{k})}^{n} = {(x_{i}, y_{i}) : 1 \leq i \leq n, i \notin F_{k}}

. Below, we first calculate the variance of

W_{1}

and then bound

ρ (W_{k_{1}}, W_{k_{2}})

with

k_{1} \neq k_{2}

. Note that

V a r (W_{k})

is equal to

\begin{matrix} \begin{matrix} V a r (\sum_{j \in F_{k}} {(f (X_{j}) + ϵ_{j} - {\hat{f}}_{(- F_{k})} (X_{j}))}^{2}) \\ = V a r (\sum_{j \in F_{k}} {(f (X_{j}) - {\hat{f}}_{(- F_{k})} (X_{j}))}^{2} + \sum_{j \in F_{k}} ϵ_{j}^{2} + 2 \sum_{j \in F_{k}} (f (X_{j}) - {\hat{f}}_{(- F_{k})} (X_{j})) ϵ_{j}) \\ = V a r (\sum_{j \in F_{k}} {(f (X_{j}) - {\hat{f}}_{(- F_{k})} (X_{j}))}^{2}) + V a r (\sum_{j \in F_{k}} ϵ_{j}^{2}) + 4 V a r (\sum_{j \in F_{k}} (f (X_{j}) - {\hat{f}}_{(- F_{k})} (X_{j})) ϵ_{j}) \\ + 2 C o v (\sum_{j \in F_{k}} {(f (X_{j}) - {\hat{f}}_{(- F_{k})} (X_{j}))}^{2}, \sum_{j \in F_{k}} ϵ_{j}^{2}) \\ + 4 C o v (\sum_{j \in F_{k}} {(f (X_{j}) - {\hat{f}}_{(- F_{k})} (X_{j}))}^{2}, \sum_{j \in F_{k}} (f (X_{j}) - {\hat{f}}_{(- F_{k})} (X_{j})) ϵ_{j}) \\ + 4 C o v (\sum_{j \in F_{k}} ϵ_{j}^{2}, \sum_{j \in F_{k}} (f (X_{j}) - {\hat{f}}_{(- F_{k})} (X_{j})) ϵ_{j}) \\ = V a r (\sum_{j \in F_{k}} {(f (X_{j}) - {\hat{f}}_{(- F_{k})} (X_{j}))}^{2}) + \sum_{j \in F_{k}} V a r (ϵ_{j}^{2}) + 4 \sum_{j \in F_{k}} V a r (ϵ_{j}) E | | f - {\hat{f}}_{(- F_{k})} {| |}_{2}^{2} \\ + 0 + 0 + 4 \sum_{j \in F_{k}} E ϵ_{j}^{3} (f (X_{j}) - {\hat{f}}_{(- F_{k})} (X_{j})) \\ = V a r (\sum_{j \in F_{k}} {(f (X_{j}) - {\hat{f}}_{(- F_{k})} (X_{j}))}^{2}) + n_{0} V a r (ϵ_{j}^{2}) + 4 n_{0} σ^{2} E | | f - {\hat{f}}_{(- F_{k})} {| |}_{2}^{2}, \end{matrix} \end{matrix}

(A3)

where we have used the assumption that the error

ϵ_{j}

is independent of

X_{j}

and

E ϵ_{j}^{3} = 0

. Now,

\begin{matrix} \begin{matrix} V a r (\sum_{j \in F_{k}} {(f (X_{j}) - {\hat{f}}_{(- F_{k})} (X_{j}))}^{2}) \\ = & E [V a r (\sum_{j \in F_{k}} {(f (X_{j}) - {\hat{f}}_{(- F_{k})} (X_{j}))}^{2} | D_{(- F_{k})}^{n})] + V a r [E (\sum_{j \in F_{k}} {(f (X_{j}) - {\hat{f}}_{(- F_{k})} (X_{j}))}^{2}) | D_{(- F_{k})}^{n}] \\ = & E n_{0} [V a r ({(f (X_{j}) - {\hat{f}}_{(- F_{k})} (X_{j}) | D_{(- F_{k})}^{n})}^{2})] + V a r (n_{0} | | f - {\hat{f}}_{(- F_{k})} {| |}_{2}^{2}) \\ = & n_{0} E (| | f - {\hat{f}}_{(- F_{k})} {| |}_{4}^{4} - | | f - {\hat{f}}_{(- F_{k})} {| |}_{2}^{4}) + n_{0}^{2} V a r (| | f - {\hat{f}}_{(- F_{k})} {| |}_{2}^{2}) . \end{matrix} \end{matrix}

(A4)

Therefore,

\begin{matrix} \begin{matrix} V a r (W_{k}) \\ = & n_{0} E (| | f - {\hat{f}}_{(- F_{k})} {| |}_{4}^{4} - | | f - {\hat{f}}_{(- F_{k})} {| |}_{2}^{4}) + n_{0}^{2} V a r (| | f - {\hat{f}}_{(- F_{k})} {| |}_{2}^{2}) \\ + & n_{0} V a r (ϵ_{j}^{2}) + 4 n_{0} σ^{2} E | | f - {\hat{f}}_{(- F_{k})} {| |}_{2}^{2} . \end{matrix} \end{matrix}

(A5)

Now for

W_{k_{1}}

and

W_{k_{2}}

k_{1} \neq k_{2}

, we calculate

C o v (W_{k_{1}}, W_{k_{2}})

. Clearly, WLOG, we can take

k_{1} = 1

k_{2} = 2

\begin{matrix} \begin{matrix} W_{1} = & \sum_{j \in F_{1}} {(f (X_{j}) - {\hat{f}}_{(- F_{1})} (X_{j}) + ϵ_{j})}^{2} \\ = & \sum_{j \in F_{1}} {(f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) + ϵ_{j} + {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) - {\hat{f}}_{(- F_{1})} (X_{j}))}^{2} \\ = & \sum_{j \in F_{1}} {(f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) + ϵ_{j})}^{2} + \sum_{j \in F_{1}} {({\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) - {\hat{f}}_{(- F_{1})} (X_{j}))}^{2} \\ + & 2 \sum_{j \in F_{1}} (f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) + ϵ_{j}) ({\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) - {\hat{f}}_{(- F_{1})} (X_{j})) \\ ≜ & C + R_{1, 1} + R_{1, 2} . \end{matrix} \end{matrix}

(A6)

Similarly,

\begin{matrix} \begin{matrix} W_{2} = & \sum_{j \in F_{2}} {(f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) + ϵ_{j})}^{2} + \sum_{j \in F_{2}} {({\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) - {\hat{f}}_{(- F_{2})} (X_{j}))}^{2} \\ + & 2 \sum_{j \in F_{2}} (f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) + ϵ_{j}) ({\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) - {\hat{f}}_{(- F_{2})} (X_{j})) \\ ≜ & D + R_{2, 1} + R_{2, 2} . \end{matrix} \end{matrix}

(A7)

Then

\begin{matrix} \begin{matrix} C o v (W_{1}, W_{2}) = & C o v (C, D) + C o v (C, R_{2, 1}) + C o v (C, R_{2, 2}) \\ + & C o v (R_{1, 1}, D) + C o v (R_{1, 1}, R_{2, 1}) + C o v (R_{1, 1}, R_{2, 2}) \\ + & C o v (R_{1, 2}, D) + C o v (R_{1, 2}, R_{2, 1}) + C o v (R_{1, 2}, R_{2, 2}) . \end{matrix} \end{matrix}

(A8)

For

C o v (C, D)

, we have

\begin{matrix} \begin{matrix} C o v (C, D) = & E C D - (E C) (E D) \\ = & E [\sum_{j \in F_{1}} {(f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) + ϵ_{j})}^{2}] [\sum_{j \in F_{2}} {(f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) + ϵ_{j})}^{2}] \\ - & {\{E \sum_{j \in F_{1}} [{(f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}))}^{2} + ϵ_{j}^{2} + 2 ϵ_{j} (f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}))]\}}^{2} . \end{matrix} \end{matrix}

For the first part, we have

\begin{matrix} \begin{matrix} E C D = & E E [C D | D_{(- F_{1}, - F_{2})}^{n}] \\ = & E E [\sum_{j \in F_{1}} [{(f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}))}^{2} + ϵ_{j}^{2} + 2 ϵ_{j} (f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}))] \\ \times & \sum_{j \in F_{2}} [{(f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}))}^{2} + ϵ_{j}^{2} + 2 ϵ_{j} (f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}))] | D_{(- F_{1}, - F_{2})}^{n}] . \\ = & n_{0}^{2} E | | f - {\hat{f}}_{(- F_{1}, - F_{2})} {| |}_{2}^{4} + n_{0}^{2} σ^{4} + 2 n_{0}^{2} σ^{2} E | | f - {\hat{f}}_{(- F_{1}, - F_{2})} {| |}_{2}^{2}, \end{matrix} \end{matrix}

and

\begin{matrix} \begin{matrix} (E C) (E D) \\ = & {(E C)}^{2} \\ = & {[n_{0} (E | | f - {\hat{f}}_{(- F_{1}, - F_{2})} {| |}_{2}^{2} + σ^{2})]}^{2} \\ = & n_{0}^{2} {(E | | f - {\hat{f}}_{(- F_{1}, - F_{2})} {| |}_{2}^{2})}^{2} + 2 n_{0}^{2} σ^{2} E | | f - {\hat{f}}_{(- F_{1}, - F_{2})} {| |}_{2}^{2} + n_{0}^{2} σ^{4} . \end{matrix} \end{matrix}

(A9)

Thus,

C o v (C, D) = n_{0}^{2} V a r (| | f - {\hat{f}}_{(- F_{1}, - F_{2})} {| |}_{2}^{2})

. For the other terms in the earlier expression of

C o v (W_{1}, W_{2})

, they can be handled similarly. Indeed,

| C o v (C, R_{2, 1}) | \leq \sqrt{V a r (C)} \sqrt{V a r (R_{2, 1})}

and the other terms are bounded similarly.

Next, we calculate

V a r (C)

. Actually, the calculation is basically the same as for

V a r (W_{k})

and we get

\begin{matrix} V a r (C) = & n_{0} E (| | f - {\hat{f}}_{(- F_{1}, - F_{2})} {| |}_{4}^{4} - | | f - {\hat{f}}_{(- F_{1}, - F_{2})} {| |}_{2}^{4}) + n_{0}^{2} V a r (| | f - {\hat{f}}_{(- F_{1}, - F_{2})} {| |}_{2}^{2}) \\ + n_{0} V a r (ϵ_{j}^{2}) + 4 n_{0} σ^{2} E | | f - {\hat{f}}_{(- F_{1}, - F_{2})} {| |}_{2}^{2} . \end{matrix}

Note that

\begin{matrix} \begin{matrix} V a r (R_{1, 1}) & = V a r (\sum_{j \in F_{1}} {({\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) - {\hat{f}}_{- F_{1}} (X_{j}))}^{2}) \\ = E (V a r (\sum_{j \in F_{1}} {({\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) - {\hat{f}}_{(- F_{1})} (X_{j}))}^{2} | D_{(- F_{1})}^{n})) \\ + V a r (E (\sum_{j \in F_{1}} {({\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) - {\hat{f}}_{(- F_{1})} (X_{j}))}^{2} | D_{(- F_{1})}^{n})) \\ = E n_{0} (| | {\hat{f}}_{(- F_{1}, - F_{2})} - {\hat{f}}_{(- F_{1})} {| |}_{4}^{4} - | | {\hat{f}}_{(- F_{1}, - F_{2})} - {\hat{f}}_{(- F_{1})} {| |}_{2}^{4}) \\ + n_{0}^{2} V a r (| | {\hat{f}}_{(- F_{1}, - F_{2})} - {\hat{f}}_{(- F_{1})} {| |}_{2}^{2}) \\ \leq n_{0} E | | {\hat{f}}_{(- F_{1}, - F_{2})} - {\hat{f}}_{(- F_{1})} {| |}_{4}^{4} + n_{0}^{2} E | | {\hat{f}}_{(- F_{1}, - F_{2})} - {\hat{f}}_{(- F_{1})} {| |}_{2}^{4} . \end{matrix} \end{matrix}

(A10)

For the other term,

\begin{matrix} \begin{matrix} V a r (R_{1, 2}) = & 4 V a r (\sum_{j \in F_{1}} (f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) + ϵ_{j}) ({\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) - {\hat{f}}_{(- F_{1})} (X_{j}))) \\ = & 4 E V a r (\sum_{j \in F_{1}} (f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) + ϵ_{j}) ({\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) - {\hat{f}}_{(- F_{1})} (X_{j})) | D_{(- F_{1})}^{n}) \\ + & 4 V a r (E (\sum_{j \in F_{1}} (f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) + ϵ_{j}) ({\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) - {\hat{f}}_{(- F_{1})} (X_{j})) | D_{(- F_{1})}^{n})) \\ = & 4 E (n_{0} V a r ((f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) + ϵ_{j}) ({\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) - {\hat{f}}_{(- F_{1})} (X_{j})) | D_{(- F_{1})}^{n})) \\ + & 4 V a r (n_{0} E (f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j})) ({\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) - {\hat{f}}_{(- F_{1})} (X_{j})) | D_{(- F_{1})}^{n}) \\ = & 4 n_{0} E V a r ((f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j})) ({\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) - {\hat{f}}_{(- F_{1})} (X_{j})) | D_{(- F_{1})}^{n}) \\ + & 4 n_{0} E σ^{2} V a r (| | {\hat{f}}_{(- F_{1}, - F_{2})} - {\hat{f}}_{(- F_{1})} {| |}_{2}^{2}) \\ + & 4 n_{0}^{2} V a r (E [(f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j})) ({\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) - {\hat{f}}_{(- F_{1})} (X_{j})) | D_{(- F_{1})}^{n}]) \\ \leq & 4 n_{0} E (E ({(f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}))}^{2} {({\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) - {\hat{f}}_{(- F_{1})} (X_{j}))}^{2} | D_{(- F_{1})}^{n})) \\ + & 4 n_{0} σ^{2} E (| | {\hat{f}}_{(- F_{1}, - F_{2})} - {\hat{f}}_{(- F_{1})} {| |}_{2}^{2}) \\ + & 4 n_{0}^{2} E (E ({(f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}))}^{2} | D_{(- F_{1})}^{n}) E ({({\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) - {\hat{f}}_{(- F_{1})} (X_{j}))}^{2} | D_{(- F_{1})}^{n})) \\ \leq & 4 n_{0} E (E {({(f (X_{j}) - {\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}))}^{4} | D_{(- F_{1})}^{n})}^{\frac{1}{2}} {(E {({\hat{f}}_{(- F_{1}, - F_{2})} (X_{j}) - {\hat{f}}_{(- F_{1})} (X_{j}))}^{4} | D_{(- F_{1})}^{n})}^{\frac{1}{2}}) \\ + & 4 n_{0} σ^{2} E | | {\hat{f}}_{(- F_{1}, - F_{2})} - {\hat{f}}_{(- F_{1})} {| |}_{2}^{2} + 4 n_{0}^{2} E (| | f - {\hat{f}}_{(- F_{1}, - F_{2})} {| |}_{2}^{2} | | {\hat{f}}_{(- F_{1}, - F_{2})} - {\hat{f}}_{(- F_{1})} {| |}_{2}^{2}) \end{matrix} \end{matrix}

(A11)

\begin{matrix} \begin{matrix} \leq & 4 n_{0} E (| | f - {\hat{f}}_{(- F_{1}, - F_{2})} {| |}_{4}^{2} | | {\hat{f}}_{(- F_{1}, - F_{2})} - {\hat{f}}_{(- F_{1})} {| |}_{4}^{2}) + 4 n_{0} σ^{2} E | | {\hat{f}}_{(- F_{1}, - F_{2})} - {\hat{f}}_{(- F_{1})} {| |}_{2}^{2} \\ + & 4 n_{0}^{2} E | | f - {\hat{f}}_{(- F_{1}, - F_{2})} {| |}_{2}^{2} | | {\hat{f}}_{(- F_{1}, - F_{2})} - {\hat{f}}_{(- F_{1})} {| |}_{2}^{2} \\ \leq & 4 n_{0} {(E | | f - {\hat{f}}_{(- F_{1}, - F_{2})} {| |}_{4}^{4})}^{\frac{1}{2}} {(E | | {\hat{f}}_{(- F_{1}, - F_{2})} - {\hat{f}}_{(- F_{1})} {| |}_{4}^{4})}^{\frac{1}{2}} + 4 n_{0} σ^{2} E | | {\hat{f}}_{(- F_{1}, - F_{2})} - {\hat{f}}_{(- F_{1})} {| |}_{2}^{2} \\ + & 4 n_{0}^{2} {(E | | f - {\hat{f}}_{(- F_{1}, - F_{2})} {| |}_{2}^{4})}^{\frac{1}{2}} {(E | | {\hat{f}}_{(- F_{1}, - F_{2})} - {\hat{f}}_{(- F_{1})} {| |}_{2}^{4})}^{\frac{1}{2}} . \end{matrix} \end{matrix}

Under the condition

E | | f - {\hat{f}}_{(δ, n)} {| |}_{4}^{4} = o (1 / n)

, it is seen that

V a r (W_{1})

and

V a r (C)

are both of order n,

V a r (R_{1, 1})

V a r (R_{1, 2})

, and

C o v (C, D)

are all of order

o (n)

, which together imply

| C o v (W_{1}, W_{2} | = o (V a r (W_{1}))

and consequently

ρ (W_{1}, W_{2}) \to 0

n \to \infty

. This completes the proof of the theorem.

References

Nan, Y.; Yang, Y. Variable selection diagnostics measures for high-dimensional regression. J. Comput. Graph. Stat. 2014, 23, 636–656. [Google Scholar] [CrossRef]
Yu, Y.; Yang, Y.; Yang, Y. Performance assessment of high-dimensional variable identification. arXiv 2017, arXiv:1704.08810. [Google Scholar]
Ye, C.; Yang, Y.; Yang, Y. Sparsity oriented importance learning for high-dimensional linear regression. J. Am. Stat. Assoc. 2018, 113, 1797–1812. [Google Scholar] [CrossRef]
Tibshirani, R. Regression shrinkage and selection via the lasso. J. R. Stat. Soc. Ser. (Methodol.) 1996, 58, 267–288. [Google Scholar] [CrossRef]
Tikhonov, A.N. On the stability of inverse problems. Dokl. Akad. Nauk SSSR 1943, 39, 195–198. [Google Scholar]
Horel, A.E. Applications of ridge analysis to regression problems. Chem. Eng. Progress. 1962, 58, 54–59. [Google Scholar]
Akaike, H. A new look at the statistical model identification. IEEE Trans. Autom. Control 1974, 19, 716–723. [Google Scholar] [CrossRef]
Hurvich, C.M.; Tsai, C.L. Regression and time series model selection in small samples. Biometrika 1989, 76, 297–307. [Google Scholar] [CrossRef]
Schwarz, G. Estimating the dimension of a model. Ann. Stat. 1978, 6, 461–464. [Google Scholar] [CrossRef]
Chen, J.; Chen, Z. Extended Bayesian information criteria for model selection with large model spaces. Biometrika 2008, 95, 759–771. [Google Scholar] [CrossRef] [Green Version]
Wang, H.; Li, R.; Tsai, C.L. Tuning parameter selectors for the smoothly clipped absolute deviation method. Biometrika 2007, 94, 553–568. [Google Scholar] [CrossRef]
Allen, D.M. The relationship between variable selection and data agumentation and a method for prediction. Technometrics 1974, 16, 125–127. [Google Scholar] [CrossRef]
Stone, M. Cross-validatory choice and assessment of statistical predictions. J. R. Stat. Soc. Ser. (Methodol.) 1974, 36, 111–133. [Google Scholar] [CrossRef]
Geisser, S. The predictive sample reuse method with applications. J. Am. Stat. Assoc. 1975, 70, 320–328. [Google Scholar] [CrossRef]
Zhang, Y.; Yang, Y. Cross-validation for selecting a model selection procedure. J. Econom. 2015, 187, 95–112. [Google Scholar] [CrossRef]
Breiman, L.; Friedman, J.H.; Olshen, R.A.; Stone, C.J. Classification and Regression Trees; Routledge: Oxfordshire, UK, 2017. [Google Scholar]
Hastie, T.; Tibshirani, R.; Friedman, J. The Elements of Statistical Learning: Data Mining, Inference, and Prediction, 2nd ed.; Springer Series in Statistics; Springer: Berlin/Heidelberg, Germany, 2009. [Google Scholar]
Yang, Y. Comparing learning methods for classification. Stat. Sin. 2006, 16, 635–657. [Google Scholar]
Harrison, D., Jr.; Rubinfeld, D.L. Hedonic housing prices and the demand for clean air. J. Environ. Econ. Manag. 1978, 5, 81–102. [Google Scholar] [CrossRef] [Green Version]
Scheetz, T.E.; Kim, K.Y.; Swiderski, R.E.; Philp, A.R.; Braun, T.A.; Knudtson, K.L.; Dorrance, A.M.; DiBona, G.F.; Huang, J.; Casavant, T.L.; et al. Regulation of gene expression in the mammalian eye and its relevance to eye disease. Proc. Natl. Acad. Sci. USA 2006, 103, 14429–14434. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Meinshausen, N.; Bühlmann, P. Stability selection. J. R. Stat. Soc. Ser. (Stat. Methodol.) 2010, 72, 417–473. [Google Scholar] [CrossRef]
Lim, C.; Yu, B. Estimation stability with cross-validation (ESCV). J. Comput. Graph. Stat. 2016, 25, 464–492. [Google Scholar] [CrossRef] [Green Version]
Yang, W.; Yang, Y. Toward an objective and reproducible model choice via variable selection deviation. Biometrics 2017, 73, 20–30. [Google Scholar] [CrossRef] [PubMed] [Green Version]

Figure 1. Probability that model with

λ_{m i n}

outperforms model with

λ_{1 s e}

for estimation.

Figure 1. Probability that model with

λ_{m i n}

outperforms model with

λ_{1 s e}

for estimation.

Figure 2. Probability that model with

λ_{m i n}

outperforms model with

λ_{1 s e}

for estimation.

Figure 2. Probability that model with

λ_{m i n}

outperforms model with

λ_{1 s e}

for estimation.