{'fraction_non_alphanumeric': 0.05046588805922522, 'fraction_numerical': 0.020550130831578275, 'mean_word_length': 4.021793125550837, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Prediction Divergence Criterion for Model Selection

The problem of model selection is inevitable in an increasingly large number of applications involving partial theoretical knowledge and vast amounts of information, like in medicine, biology or economics. The associated techniques are intended to determine which variables are"important"to"explain a phenomenon under investigation. The terms"important"and"explain"can have very different meanings according to the context and, in fact, model selection can be applied to any situation where one tries to balance variability with complexity. In this paper, we introduce a new class of error measures and of model selection criteria, to which many well know selection criteria belong. Moreover, this class enables us to derive a novel criterion, based on a divergence measure between the predictions produced by two nested models, called the Prediction Divergence Criterion (PDC). Our selection procedure is developed for linear regression models, but has the potential to be extended to other models. We demonstrate that, under some regularity conditions, it is asymptotically loss efficient and can also be consistent. In the linear case, the PDC is a counterpart to Mallow's Cp but with a lower asymptotic probability of overfitting. In a case study and by means of simulations, the PDC is shown to be particularly well suited in"sparse"settings with correlated covariates which we believe to be common in real applications.


Introduction
Model selection is an important and challenging problem in statistics. Indeed, it becomes unavoidable in more and more applications involving incomplete theoretical knowledge about the phenomenon under investigation and important amounts of available information, like in medicine, biology, economics, etc. Very often model selection is about choosing among a set of predictors, the subset that best predicts or explains a response variable.
In general, in regression problems, the task of model selection is performed using either stepwise approaches or through penalized estimation based methods. The former optimizes a criterion such as for example the AIC on a given sequence of models. The latter produces simultaneously a sequence on which a criterion is employed to select a model, the lasso being a popular example. More details on model selection methods are provided a literature review given in Appendix A.
One of the main challenges in model selection is the development of simple and rapid techniques that can deal with large data sets especially in sparse settings, i.e. when only a small minority of the variables at hand are significant predictors. This is for example the case in areas such as genetics or economics and the available information can now easily contain millions of observations and hundreds of thousands of potential predictors. Even moderately large dataset can lead to very large ones when considering for instance all potential interactions. In this context, stepwise techniques can often provide improvements in term of computational speed and/or prediction accuracy (see e.g. Foster and Stine, 2004;Lin et al., 2011). Moreover, a stepwise approach is also often considered as more suitable when there is a natural sequence of potential predictors. Such a sequence can arise from an priori theoretical insight and is the one of interest for the scientist, like in experimental settings where main effects are more important than interactions (for mixed models, see e.g. Müller et al., 2013). Another example comes from natural sciences were a sequence of models corresponds to a sequence of null hypotheses of interest that need to be tested and upon which a model selection approach is used (see e.g. Johnson and Omland, 2004).
In this paper, we propose a new class of model selection criteria that are suited for stepwise approaches or can be used as selection criteria in penalized estimation based methods. This new class, called the d-class of error measure, generalizes Efron's q-class (see Efron, 1986). This class not only contains classical criteria such as Mallow's C p or the AIC, but also enables one to define new criteria that are more general. Within this new class, we propose a model selection criterion based on a prediction divergence between two nested models' predictions that we call the Prediction Divergence Criterion (PDC). The PDC provides a different measure of prediction error than a criterion, say C, associated to each potential model within a sequence and for which the selection decision is based on the sign of ∆C. The PDC directly measures the prediction error divergence between two nested models and provides different criteria than ∆C. As an example, we consider linear regression models and propose a PDC criterion that is the direct counterpart of Mallow's C p . We show that a selection procedure based on the PDC, compared to the C p , has a smaller probability of overfitting and a negligible asymptotic probability of selecting a larger model for models with more than one additional non significant covariate (which is not the case for selection procedures based on the C p ).
In practice, the PCD can be applied to a known or estimated sequence of models to select candidate models. Compared to other criteria used in stepwise approaches the PDC appears to often perform better in finite samples, especially in sparse settings. Both in simulations and a real example we find that the PDC tends to select far smaller models with generally better or comparable out-of-sample predictive performances. Compared to penalized estimation based methods, the PDC, even based on a rather simple ordering rule, appears to selects at least as many of the significant variables but with far fewer of the non-significant ones hence providing smaller models. In terms of prediction accuracy the PDC seems to achieve comparable results to adaptive penalized estimation based methods in sparse settings.
This point is illustrated through the analysis of a dataset on childhood malnutrition in Zambia in Section 2. We compare the performance of several model selection procedures and conclude that different procedures can lead to very different selected models, which in practice rises the problem of choosing which model selection procedure to use. The simulation study in Section 5, which compares the finite sample performance of the PDC's estimator with other model selection techniques (including adaptive ones), reveals that in sparse and correlated covariates settings, the PDC tends to select as many as the significant variables than the other best performing methods, but far less of the noisy ones, with comparable out-of-sample prediction errors. This leads then to conclude, for the analyzed dataset, that the PDC is the most reliable selection criterion that provides the smallest model including most of the explanatory variables.
The reminder of the paper is organized as follows. In Section 3 we introduce the d-class of error measures which generalizes Efron's q-class and we derive the associated optimism theorem. We also present a new class of model selection criteria, which is based on the d-class of error measures (i.e. the PDC), and derive a corresponding unbiased estimator. In Appendix C, for the linear model and L 2 divergence, we compare the PDC to Mallows' C p and also show that the latter is a particular case of the the PDC class. In Section 4, we consider the PDC for variable selection in linear regression models and derive the asymptotic properties, namely the probabilities of overfitting, consistency, asymptotic loss efficiency and the Signal to Noise Ration (SNR) (see e.g. McQuarrie and Tsai, 1998). We also propose, for these models, a suitable sequence of nested models that we use in our simulation exercise. The use of the PDC to other models such as smoothing splines or for the order of autoregressive models can be found in Guerrier (2013).

An Illustrative Example: Childhood Malnutrition in Zambia
Childhood malnutrition is considered to be one of the worst health problems in developing countries (e.g. United Nations Children 's Fund, 1998). Both a manifestation and a cause of poverty, malnutrition is thought to contribute to over a third of death in children under five years old globally (United Nations Children's Fund, 2012). Moreover, it is well established in the medical literature that maternal and child undernutrition have considerable consequences for adult health and human capital (see e.g. Victora et al., 2008 and the references therein). Such conditions are, for example, associated with less schooling, reduced economic productivity, and for women lower offspring birthweight. It has also been reported that lower birthweight and undernutrition in childhood have an influence on cancer occurrence and are risk factors for high glucose concentrations, blood pressure, and harmful lipid profiles.
The evolution of (childhood) malnutrition is monitored through regular demographic and health surveys conducted by Macro International in cooperation with the World Health Organization. These data have been made publicly available and can be obtained from www.measuredhs.com. We will consider here a dataset obtained from a demographic and health survey conducted in Zambia in 2007 (see Zambia DHS, 2007 for details).
Undernutrition is generally assessed by comparing anthropometric indicators such as height or weight at a certain age to a reference population. A well established measurement for the study of acute malnutrition is given by (see World Health Organisation, 1995 for details): where H i,j , µ j and σ j denote, respectively, the height of the i th child at age j, the median height of a child of the same age in the reference population and the associated standard deviation. Several variables are assumed to have a determinant influence on undernutrition. These are generally linked to education, income, and nutritional situation of the parents, access to clean water and sanitation, and primary health care, and immunization facilities (see e.g. Madise et al., 1999). Therefore, taking into account the available information in the considered dataset, we selected the covariates given in Appendix B to explain Y i .
After removing missing data, the dataset contains n = 1927 observations and p = 31 covariates. In a first step, we compare the performance of the lasso, the adaptive lasso, MCP, SCAD, and the forward stepwise AIC, BIC (see Appendix A for more detailed explanations) and PDC (based on the L 2 divergence, with λ n = 2 and using the ordering rule given in (25), see Sections 3 and 4). The selection criteria as well as the specifications for the lasso and adaptive methods are presented in Appendix F. The selected models are compared using the mean squared and median absolute error of prediction (MSPE and MAPE), as measured by 10-fold cross-validation. That is, we split the data into ten roughly equal-sized parts and for each part, we carry out model selection using the other nine parts of the data and calculate the MSPE and MAPE of the chosen model when predicting the k th part of the data, i.e.
where || · || p denotes the L p norm, y (k) the k th split of the vector of observation y andŷ (k) the prediction of y (k) made without y (k) . For all methods, the data were split in the same way.
In a second analysis, we include in the model all possible first order interactions between the regressors, hence leading to p = 496 potential regressors.
In Table 1 the MSPE and MAPE for the different selection methods are presented. One first notices that the performances, as measured by the MSPE or by the MAPE, are the same across the different selection procedures, for both performance measures. We also see that the PDC is the selection procedure producing overall and by far the smallest model. Regularization methods need to included many more variables in the model to achieve the same prediction accuracy as the PDC and other stepwise procedures considered here. With the addition of interactions as potential explanatory variables, two adaptive methods fail in providing a selected model (MCP and SCAD).
However, the analysis of a real example does not necessarily provide the whole story. In particular, it is important for the model selection procedures not to miss important information (i.e. significant variables). In the simulation study provided in Section 5, we find that with correlated covariates and in sparse settings, the PDC (with a simple ordering rule), compared to other methods, not only chooses at least comparable proportions of the significant variables, but also selects less of the noisy ones.

The d-Class Error Measures and an Associated Prediction Divergence Criterion
In this section, we develop a general framework that provides model selection criterion estimators in a rather straightforward manner. This framework generalizes Efron's q-class of error measures (see Efron, 1986).
Consider a random variable Y distributed according to model F θ , possibly conditionally on a set of fixed covariates x = [x 1 . . . x p ]. We observe a random sample Y = (Y i ) i=1,...,n supposedly generated from F θ , possibly together with a non-random n × K, K ≥ p matrix of inputs X. Given a prediction functionŶ that depends on the chosen model, Efron (1986) uses a function Q(·, ·) based on the q-class error measure to build an out-of-sample "performance" measure of a given prediction rule. The q-class of error measures is given by whereq(v) is the derivative of q(·) evaluated at v. The particular choice of q(u) = u(1 − u) gives the squared loss function Q(u, v) = (u − v) 2 . The prediction error measure is quantified by the (out-of-sample) expected prediction error ..,n a random variable distributed as Y , which can be interpreted as an out-of-sample version of Y . As throughout this article, E[·] and E 0 [·], denote expectations under the distribution of Y i |x i , respectively Y 0 i |x i , the correct model. The expectations are, depending on the context, simple or multiple.
Efron's optimism theorem (see Efron, 2004) demonstrates that Hence, an estimator of EPErr is obtained as where, depending on the distribution of Y i |x i , cov(·, ·) is obtained analytically up to a value of θ, the model's parameters, which is then replaced byθ, or by resampling methods (see e.g. Efron, 2004).
We may extend this methodology and the q-class of error measures by changing Q(Y 0 i ,Ŷ i ) in (3) to a more general class of error measures D(·, ·) between two equidimentional vector valued functions g 1 Y 0 ,θ 1 and g 2 Y ,θ 2 whereθ 1 andθ 2 denote the estimated parameter vectors associated, respectively, to the models F θ1 and F θ2 . Such a criterion can be defined without loss of generality as where the expectation is multidimensional. The divergence D(g 1 (Y 0 ,θ 1 ), g 2 (Y ,θ 2 )) is said to belong to the d-class of error measures if the function D(·, ·) is a valid Bregman divergence and if the functions g j (θ j , Y ) are equidimensional and associated, respectively, to the models F θj , j = 1, 2. In some sense, Bregman divergences are the multivariate equivalent of Efron's q-class (see Bregman, 1967 for more details). This divergence encompasses squared error, relative entropy, logistic loss, Mahalanobis distance and other error measures. It has often been used as loss function in the context of model selection (see e.g. Zhang, 2008;Zhang et al., 2009Zhang et al., , 2010. The Bregman divergence between two equidimensional vectors u and v (of the same dimension) is defined as where ψ(·) is a scalar and ∇ψ(v) represents the gradient vector of ψ(·) evaluated at v. The function ψ(·) is strictly convex and differentiable. For example, a squared loss function D(u, v) = ||u − v|| 2 2 is obtained when Based on the d-class of error measure and for the criterion defined in (5) we can derive the following "optimism" theorem whose proof is presented in Appendix D.1.
Theorem 1: Let the divergence D (·, ·) be a valid d-class error measure based on ψ (·) and assume that Then The direct consequence of Theorem 1 is that for any criterion C as defined in (5) one can construct an estimator, say C similarly to (4) for the q-class of error measures, using (3). Indeed a "natural" (and unbiased) estimator of C is C = D g 1 y,θ 1 , g 2 y,θ 2 + tr cov g 1 y,θ 1 , ∇ψ g 2 y,θ 2 where as in (4), depending on the distribution of Y |X, cov(·, ·) is obtained analytically up to a value of θ 1 and θ 2 , the models' parameters, which are then replaced byθ 1 andθ 2 , or by resampling methods (see e.g. Efron, 2004). We shall refer to the first and second terms of (7) as the apparent divergence and the divergence optimism, respectively.
The class (7) encompasses many selection criteria such as Mallows' C p which is presented as an example in Appendix C together with meaningful properties a criterion should satisfy for the task of model selection.
As a model selection criterion to choose between two models, say M j nested in M k , we propose to consider the class C in (5) with g 1 Y 0 ,θ j =Ŷ 0 j and g 1 Y ,θ k =Ŷ k . This class directly compares the out-of-sample prediction computed in the smaller modelŶ 0 j with the in-sample prediction in the larger modelŶ k , quantified by the Bregman divergence D (·, ·) (based on ψ(·)).
More formally, a criterion that compares the prediction divergence between two models belongs to the following class: where the expectation is multidimensional. By using Theorem 1 as in (7) we obtain D (ŷ j ,ŷ k ) + tr { cov [ŷ j , ∇ψ (ŷ k )]} as an unbiased estimator of the PDC j,k . Although this estimator can be computed analytically or using resampling methods for any valid Bregman divergence, we shall only consider here the squared loss function (the same as the one for the C p ) to get In some sense, PDC j,k is comparable (but not equal) to Mallow's C p in the PDC class since both criteria are based on the same loss function. It will thus be of particular interest to compare (9) with the C p to understand the differences between the PDC approach and the classical model selection approach. This comparison is presented in Section 4 in the context of the linear regression model.
The PDC has an intuitive interpretation. Indeed, if the smaller model is not correct, the additional elements in the larger model create differences in the predictions and therefore should be accounted for. For a "suitable" sequence of nested models with increasing complexity and such that model M j is nested in model M j+1 , PDC j,j+1 is expected to be minimal when j =  0 < K, K being the number of potential sequentially nested models, and M 0 denotes the correct (or closest to the correct) model. Indeed, while j <  0 , we expect PDC j,j+1 to be relatively large since model M j is missing some elements of the correct model M 0 which are included in model M j+1 . This is also true with PDC j,j+m , j < j + m ≤ K or PDC j−m,j , 1 < j − m < j. On the other hand, if j ≥  0 , PDC j,j+1 (or indeed PDC j,j+m , m > 0) is relatively small compared to when j <  0 because both models include the correct one. Among all models j ≥  0 , PDC j,j+m should be minimised at j =  0 and m = 1 since PDC 0,0+1 compares the prediction of the correct model with the least overfitted one.
In the case of the linear regression model, we derive in Section 4 the (asymptotic) properties of the PDC j,j+1 . In particular, we show in Theorem 2 that for sufficiently large sample size n we have that E PDC 0,0+1 ≤ E PDC j,j+m for j and m such that 0 < j < K, m > 0 and j + m ≤ K + 1. We also have that E PDC 0,0+1 = E PDC j,j+m if and only if j =  0 and m = 1. This confirms the intuitive explanation given above.
Many authors (see e.g. Bhansali and Downham, 1977) have examined the penalty function of the AIC (and of other criteria) and defined, for example, the AIC α in which the term 2 (see Table A-5 in Appendix F) of the conventional AIC is replaced by α. We follow this strategy and define PDC where λ n is a constant depending possibly on the sample size n.
Hence, assuming that there exist K competing nested models (and that the largest model is not the correct one) for describing the behavior of Y , we propose to choose the model M λn satisfying If a clear sequence of competing nested models does not exist, one can build one prior to applying the selection rule (11). This will be explained when treating the linear regression model in Section 4.4 (in particular see iterative rule (25)).

Stepwise Forward Selection for the Linear Regression Model
We consider in this section the model given in the following setting.
Setting A: Consider a random variable Y distributed according to a model that generates random samples y = (y i ) i=1,...,n through In the linear regression setting, we also observe a non-random n × K full rank matrix of inputs X, such that lim n→∞ 1/nX T X = M with M a positive definite matrix. We try to approximate (12) by means of the linear model with θ = β ∈ B ⊆ K . We actually consider in this setting the case where the covariates are "ordered". We hence rewrite β := β j = [β 1 , ..., β j , 0, ..., 0] and define β 0 = [0, ..., 0]. The columns of X are ordered accordingly.
The Least Squares Estimator (LSE) of β is given bŷ For each possible model, say M o constructed from (13), we associate (sample) predictionŝ The PDC λn j,k defined in (10) can be simplified for two linear nested candidate models. Indeed, let model M j be nested within model M k and so that dim(β j ) = j < dim(β k ) = k. Then, the PDC λn j,k based on the squared loss function defined by ψ(z) = z T z, for λ n = 2 (then PDC 2 j,k ≡ PDC j,k ), is equal to Note that (15) is not equal (or proportional) to the difference in Mallows' C p computed at respectively models M j and M k .
We obtain the modified version PDC λn j,k as: When the value of σ 2 ε is unknown, one may replace σ 2 ε by a consistent estimator, sayσ 2 ε , like the LSE at the full (possible) model. In Theorem 2 below, we show that PDC λn j,j+m is expected, for sufficiently large sample size, to reach its smallest value for j =  0 and m = 1, where  0 is the number of covariates in the correct underlying model M 0 . This motivates the selection rule defined in (11). The proof of Theorem 2 is presented in Appendix D.2.
Theorem 2: We assume Setting A and that σ 2 ε is either known or estimated by a consistent estimator, sayσ 2 ε . Then, for j and m such that 0 < j < K, m > 0 and j + m ≤ K + 1 we have that for sufficiently large n We also have that E PDC 0,0+1 = E PDC j,j+m if and only if j =  0 and m = 1.
Given Setting A, the selection problem hence amounts at selecting the number j. Let the class of models J correspond to all possible models that can be constructed. Then we define the "best" model  0 ∈ J as whereβ j denotes the j th element ofβ K . We also define the class of models J 0 which "includes"  0 and we say that model j ∈ J 0 if j ≥  0 .
Remark A: The definition of  0 in (18) may not always be a suitable definition of the model we should select. Indeed, in some situations it could be possible that some elements of the vector β are decreasing functions of n which vanish as n → ∞. Therefore, in such a setting, (18) should not be employed and a more complex definition of  0 should be used. However, in this article we will assume that none of the elements of β depend on the sample size n.
Suppose that  0 is estimated using (16) together with the selection rule given in (11), we obtain To derive the asymptotic property of the selection procedure given in (19), we will consider the following assumptions: (A.1) When σ 2 ε is unknown it can be replaced by a consistent estimatorσ 2 ε . (A.2) The scalars λ n and  0 are such that λ n  0 = O ( √ n) and λ n > 0.
Remark B: Assumption (A.1) will be employed in nearly all the results considered in this article. In practice, a popular choice ofσ 2 ε is given byσ 2 (i.e. the LSE at the full model) which is consistent under model (12) with µ = Xβ 0 . However, this estimator is not consistent if µ = Xβ 0 . An alternative could for instance be to use of the nearest neighbor method (see e.g. Stone, 1977). A more detailed discussion on this matter can be found in Section 4 of Shao (1997).
Remark C: Assumption (A.2) imposes some conditions on λ n (as presented in (19)) and on  0 . In most situations, it is quite reasonable to assume that  0 does not depend on the sample size n. However, this is not always the case (see e.g. the "one-mean versus k-mean" example of Shao (1997)). The penalty term λ n can also be an increasing function of n and we will see in Theorem 5 that this is a required condition for the selection procedure (19) to be consistent. However the product of λ n and  0 can never increase at a faster rate than √ n.
In the following subsections we study in turn the probabilities of underfitting and overfitting as well as the conditions for consistency and asymptotic loss efficiency for PDC λn j,k . We also study the SNR of PDC λn j,k and compare it to other model section criteria in Theorem A-8 given in Appendix E.

Underfitting
The least we can expect from any reliable model selection procedure is that for sufficiently large n the probability of selecting a model that does not belong to J 0 (i.e. the probability of underfitting) tends to zero. Under reasonable conditions (see e.g. Remark A) almost all model selection criteria such as the AIC, Mallow's C p or the BIC have a nil asymptotic probability of underfitting. Theorem 3 examines the asymptotic probability of selecting a model that does not belong to J 0 using λn as defined in (19). The proof of Theorem 3 is presented in Appendix D.3.

Probability of Overfitting
The next theorem (whose proof is presented in Appendix D.4) determines the asymptotic probability of overfitting in the case where we only consider choosing between models  0 and j ∈ J 0 \ { 0 }. It is based on the variance-gamma distribution (see e.g. Kotz et al., 2001) which is described in Appendix G. In this situation, model  0 will be selected if the following inequality is true: where m = j −  0 and j ∈ J 0 \ { 0 }. In general, we only consider the case where m = 1 as in the selection rule given in (19), but we provide here the results for any m ∈ N + such that j + m ≤ K + 1.
Theorem 4: Let j ∈ J 0 \ { 0 } and let G m be a random variable following a modified variancegamma distribution with parameter m = j −  0 . We assume that j + m ≤ K + 1 and under Setting A and Assumption (A.1) we have that Remark D: In the case where σ 2 ε is known the result of Theorem 4 is no longer asymptotic and the limit can therefore be removed. This directly follows from (A-25) in Appendix D.4.
Theorem 4 enables to compare the asymptotic probabilities of overfitting by m variables of PDC λn j,j+1 with other approaches such as Mallow's C p . The results are presented in Table 2 and are based on the derivation presented in McQuarrie and Tsai (1998, Section 2.5). It can be observed that PDC j,j+1 (i.e. λ n = 2 resulting from Theorem 1) will overfit less than the C p (and other asymptotically equivalent methods such as the AIC or the FPE) but more than the AICu. Clearly, PDC λn j,j+1 based on a λ n that tends to infinity with n has, similarly to the BIC, HQ and HQc, nil asymptotic probabilities of overfitting.
In Appendix D.5 we derive the asymptotic probability of overfitting for PDC λn j,j+1 (see Theorem A-7). Woodroofe (1982) and Zhang (1992) showed that the C p (and asymptotically equivalent methods) are such that For PDC j,j+1 , the counterpart of Mallow's C p in the PDC class , we have that The PDC j,j+1 has therefore a larger probability of selecting the "best" model compared to the C p or the AIC. Figure 1 presents the limiting overfitting probabilities for Mallow's C p and PDC j,j+1 . It can be observed that unlike the C p , PDC j,j+1 is not much affected by the number of models larger than  0 (i.e. K −  0 ) and tends to select less overfitted models. Indeed, when there is only one model larger than  0 , PDC j,j+1 adds on average about 0.10 additional variables in the model, and when the number of models tends to infinity, the number of additional variables only increases to about 0.11. Moreover, regardless of the number of candidate models, the 3.52 % 0.14 % 0.32 % 10 2.93 % 0.09 % 0.21 % asymptotic probability for PDC j,j+1 to select a model with three or more additional variables is virtually nil. From Corollary A-3 in Appendix D.6 (see Remark G), we show that PDC j,j+1 selects on average 0.11 more variables, while for the C p (and the AIC) this number is 0.946 (see also Woodroofe, 1982 andZhang, 1992 for details).

Consistency and Asymptotic Loss Efficiency
A selection procedure which uses to estimate  0 is said to be consistent if Note that (22) implies Pr (L n () = L n ( 0 )) n the squared error loss for a given M j producing a response predictionŷ j and µ is the vector of When the true model is finite, the BIC is consistent (see e.g. Haughton, 1988), while the AIC (and the C p which are asymptotically equivalent as shown e.g. in Nishii, 1984) have a non-nil (asymptotic) probability of overfitting. In finite samples, consistency is not necessarily a good property since consistent selection criteria may tend to underfit, so that the chosen models could have larger prediction errors.
In some cases, a selection procedure does not fulfil (22) but is still "close" to  0 in the following sense that is weaker than (22): A selection procedure satisfying (23) is said to be asymptotically loss efficient (see Shao, 1997).
The next theorem enables to relate PDC λn j,j+1 to properties (22) and (23) which are without surprise linked to the choice of λ n . The proof of this theorem is presented in Appendix D.6.
Theorem 5: Under Setting A and Assumptions (A.1) and (A.2), we have that If in addition λ n is such that lim n→∞ λ n = ∞ then we also have that

Ordering rule
In order to apply rule (19), one needs a sequence of competing nested models that includes model  0 . In practice, this sequence might not always exist naturally, and we have instead a set of K predictors which can generate a very large number of nested sequences. We propose here an algorithm for finding such a sequence, and will show in Theorem 6 (below), that this algorithm orders the variable in a suitable sequence for sufficiently large n.   (Woodroofe, 1982,  We suppose that a "null" model M 0 of say size 0 ≤ r < K which represents the smallest possible model nested in the correct one is available. Such a model is typically the model Y = ε (in that case, r = 0). From M 0 , we choose model M 1 that includes a single additional predictor chosen among the K − r available ones, such that the PDC λn between the null model and the model with the chosen predictor is maximal. Then, starting from M 1 , one applies again the same rule until a sequence of K − r nested models is obtained. This allows one to have a convex optimization function for the PDC λn .
Let PDC λn (k) j,j+1 denote the PDC estimator between M j and M j+1 with the latter containing the additional predictor k among the K − j available ones, then at model M j , the rule is It can be noted that for the squared loss function, maximizing PDC (19) is equivalent to simply maximizing the apparent divergence since all models have the same divergence optimism (i.e. λ n σ 2 ε j). Therefore, the rule given in (24) Using (25) instead of (24) has the advantage of allowing to "order" the covariates without an estimate of σ 2 ε .
Moreover, rule (25) actually provides the same order as the one based on the C p , i.e the residual sum of squares. Indeed, we have that In Theorem 6 (below) we verify that (25) enables to correctly order the covariates. The proof of this result is given in Appendix D.7.
Theorem 6: In Setting A, let X denote the matrix X whose columns are reorganised according to the iterative rule (25). Then the first  0 columns of X contain all significant elements of β in the sense that (see (18)) Remark E: In finite samples, the iterative rule (25) together with the PDC λn j,j+1 selection criterion defined in (19) is likely to produce "sparse" models. Indeed, the rule (25) maximizes at each step the apparent divergence while (19) adds at each step a larger penalty (since the number of parameters increases). Thus, the apparent divergence is expected to decrease at each iteration while the optimism penalty increases. In finite samples, there certainly exists, in some situations, another permutation of the columns of the design matrix X, say X • (which is different from X ) that leads to a model with more covariates than the model that would be obtained by applying (25) and then (19).
Remark F: Other ordering rules can in principle be applied to obtain an ordered design matrix X whose first  0 columns contain (asymptotically) all significant elements of β. The lasso sequence is such a possibility (see Donoho and Elad, 2002;Donoho and Huo, 2002;Donoho, 2006 andZou, 2006). The comparison with our procedure is however left for future research.

Simulation Study
In this section we consider three different simulation settings that corresponds to three different possible situations encountered in practice. The aim is to compare the behavior of different selection rules to, on the one hand, study how the asymptotic properties apply in finite samples, and, on the other hand, to study the effect of different settings on the performance of the different selection rules.
We consider the PDC λn j,j+1 (defined in (19)) using three different λ n , namely λ n = 2 ( PDC j,j+1 ), λ n = log (n) ( PDC • j,j+1 ) and λ n = 2 log (log (n)) ( PDC j,j+1 ). We compare PDC λn j,j+1 with other selection methods in the linear regression framework. The latter include the stepwise forward approach using explicit model selection criteria such as the AIC, AICc, AICu, FPE, FPEu, BIC, HQ, HQc (see Table A-5 in Appendix F), which are applied to a prior ordering of the variables. This ordering processes use rule (24) in which PDC λn j,j+1 is replaced by the corresponding criterion (i.e. AIC, AICc, etc.). We also include the lasso which provides an ordered list of variables and the model is chosen by means of the C p statistic, as well more sophisticated penalized estimation methods such as the elastic net (enet), the adaptive lasso (a-lasso), the MCP and the SCAD. We also compute the LSE on the complete model as a benchmark. The packages used and the chosen settings for these methods, as well as for the lasso, are provided in Appendix F. Adaptive methods are expected to perform better because they make extensive use of the model assumptions. It should be noted that the PDC procedure could also be extended to an adaptive one, but this is left for further research.
The performance is measured by means of the indicators provided in Table 3 as well as distributions (boxplot) of the mean squared error of prediction (estimated on a test set), i.e and the estimation precision assessed using the average mean squared error of estimation (estimated in the training set only), i.e The boxplots and the medians of both PE y and MSE β are then used to compare the methods. Across simulation settings, we vary the the level of sparsity, the level of correlation among the covariates as well as the SNR of the slopes.
Simulation 5.1. In this simulation we consider a sparse setting where we expect the PDC λn j,j+1 criteria to perform well. We consider a linear model with  (26)) computed on test samples. (27)) computed on test samples. and σ 2 = 1. The covariates are standard normal realizations with pairwise correlations between x j and x k (i.e. j th and k th columns of X) arbitrarily set to corr(x j , x k ) = 0.5 |j−k| . This situation corresponds to a theoretical R 2 of 66.7% and to a SNR for the significant slope parameter of 2.0. We chose n = 80 for the training and n = 800 for the test samples.
Simulation 5.2. This simulation considers a relatively "dense" setting where we expect the PDC λn j,j+1 criteria to perform poorly. We consider a linear model with and σ 2 = 1. The covariates are standard normal realizations with pairwise correlations between x j and x k arbitrarily set to corr(x j , x k ) = 0.75 |j−k| . This situation corresponds to a theoretical R 2 of 52.6% and to a SNR for the significant slope parameter of 1.1. We used n = 100 for the training and n = 10 3 for the test samples.
Simulation 5.3. In this simulation we consider a relatively sparse setting but where some of the elements of β are small and thus difficult to identify. We expect the PDC λn j,j+1 criteria to perform well in terms of prediction and estimation error but not necessarily well in terms of model identification. We consider a linear model with β = (2, 0, 1, 2, 0, 1, 0, ..., 0 16 , 0.1, ..., 0.1 6 , 0, ..., 0 16 , 2, 0, 1, 2, 0, 1) and σ 2 = 4. The covariates are standard normal realizations with pairwise correlations between x j and x k arbitrarily set to corr(x j , x k ) = 0.5 |j−k| . This situation corresponds to a theoretical R 2 of 88.5% and to a SNR of approximately 7.7. We choose n = 100 for the training and n = 10 3 for the test samples. Table 4 presents the performance of some model selection criteria using the measures in Table 3 for respectively the three different settings. The complete results are given in Tables A-6 to A-8 in Appendix I. Figures 2 to 4 show the performance of all criteria in terms of prediction and estimation error.
When considering other performance criteria in Table 4 (see also Tables A-6 to A-8 in Appendix I), in the very sparse setting (Simulation 5.1) the proportion of times the PDC selects the correct model is of 56%, and the adaptive lasso and MCP have a slightly better performance. The elastic net, SCAD and lasso are significantly worse, and the other selection criteria are just way out. By considering the proportion of times the correct model is nested within the selected model (Inc. [%]), one can see that at least all methods always include the correct model, with the notable exception of the elastic net. The number of extra noisy variables that are included in the selected model (false+) show that the PDC, adaptive lasso and also the MCP are the most accurate model selection procedures in catching the correct model when the setting is very sparse. In the dense setting (Simulation 5.2), the proportion of times the correct model is selected is very low overall. The enet and adaptive methods include the true model less than half the times, but perform better than all other methods. In this setting, we can say that none of the methods have a really satisfactory performance, and the LSE at the full model seems to be sufficient. Finally, in the low signal setting (Simulation 5.3), without surprise, no method selects the correct model most of the times. The average number of true positives (significant variables) included in the selected model is comparable across methods, while the average number of false negative (non significant variables) included in the model is the lowest for the PDC.
We hence can conclude that in sparse settings, the PDC, with a simple ordering rule, tends to have a similar performance not only in terms of out-of-sample prediction error but also in choosing the correct variables when compared to the other best methods, while it includes the smallest number of noisy variables. This might suggest that with a real dataset, the PDC will pick up at least as many of the significant variables as the other methods, but will pick up only a few of the noisy ones. This is a probable explanation of the conclusions drawn in the data analysis of the Childhood Malnutrition in Zambia in Section 2. Moreover, the PDC is simple to implement and hence can easily be used on very large datasets, without prior screening or subset selection.  (26) and (27) (26) and (27) (26) and (27)

Appendices A Model selection methods
A common model selection procedure consists in computing a criterion, say C, associated to a suitable sequence of potential models, and in choosing the one(s) that optimize this criterion. Many criteria have been proposed and the most popular ones include Mallow's C p (Mallows, 1973) based on prediction error, Akaike's Information Criterion (AIC) (Akaike, 1974), based on the Kullback-Leibler divergence, and the Bayesian Information Criterion (BIC) (Schwarz, 1978).
Beside the C p , the AIC and the BIC a great number of criteria have been proposed in the literature (see e.g. Hannan and Quinn, 1979;Foster and George, 1994;Zheng and Loh, 1997;Tibshirani and Knight, 1999;George (2000) and the references therein). Kashid and Kulkarni (2002) consider a criterion based on the difference in prediction error. Alternative formulas for the AIC penalty have been proposed (see e.g. Hurvich and Tsai, 1989;McQuarrie et al., 1997;McQuarrie and Tsai, 1998 and Table A-5 in Appendix F). An alternative to explicit model selection criteria is to estimate the expected out-of-sample prediction error by simulation methods such as bootstrap (see e.g. Efron, 1983) or cross validation (see e.g. Shao, 1993;Reiss et al., 2012). These methods are strongly linked to explicit model selection criteria; for example, Shao (1993) showed that the leave-one-out cross validation is asymptotically equivalent to the AIC. Moreover, Efron (2004) demonstrated that model-based selection criteria such as the AIC have better model selection performance compared to nonparametric methods like cross validation, assuming that the model is believable.
Other methods that are nowadays well established are based on penalized estimation and include the bridge penalty type (Frank and Friedman, 1993), the nonnegative garrote (Breiman, 1995), the Smoothly Clipped Absolute Deviation or SCAD (Fan and Li, 2001), the elastic-net (Zou and Hastie, 2005), the Dantzig Selector (Candès and Tao, 2007), a correlation based penalty (Tutz and Ulbricht, 2009) or the Minimax Concave Penalty or MCP (Zhang, 2010). Some recent work propose extensions of these methods to improve in some sense their performance. Zou (2006) propose an adaptive lasso, where adaptive weights are used for penalizing different coefficients, so that the resulting model selection method enjoys the oracle property. Meinshausen (2007) proposes the relaxed lasso, a two-stage procedure that provides a continuum of solutions that include both soft-and hard-thresholding of estimators. Zou and Zhang (2009) propose the adaptive elastic-net, that combines the strengths of the quadratic regularization and the adaptively weighted lasso shrinkage. Meinshausen and Bühlmann (2010) uses a very generic subsampling approach to determine the amount of regularization in the lasso such that a certain familywise type I error rate in multiple testing can be conservatively controlled for finite sample size. Sun and Zhang (2012), based on Staädler et al. (2010), propose to jointly estimate the slope and the noise level in the regression model using a penalized estimating equation (scaled lasso). Adaptive methods have also been proposed for computing the penalty of explicit selection criteria. Shen and Ye (2002) propose an adaptive model selection procedure that uses a data-adaptive complexity penalty based on a concept of generalized degrees of freedom. Müller and Welsh (2010) consider penalty curves, in which the penalty multiplier is varied over a given interval, for improving the performance of loss based criteria such as the AIC.
Finally, model selection tools for the regression model that are robust to small departures from the model (e.g. outlying observations) have also been proposed as alternatives to either models selection criteria or to robustify penalized regression method. Ronchetti (1982) propose a robust AIC, Machado (1993) a robust BIC, Ronchetti (1997)

B Childhood Malnutrition in Zambia
The selected explanatory variables to explain the response score (1) in the Child Malnutrition in Zambia dataset, taken from www.measuredhs.com, are the following: • breastfeeding duration (month); • age of the child (month); • age of the mother (years); • Body Mass Index (BMI) of the mother (kg/meter 2 ); • height of the mother (meter); • weight of the mother (kg); • region of residence (9 levels: Central, Copperbelt, Eastern, Luapula, Lusaka, Northern, Northwestern, Southern and Western); • mother's highest education level attended (4 levels: No education, Primary, Secondary and Higher); • wealth index factor score; • weight of child at birth (kg) ; • sex of the child; • interval between the current birth and the previous birth (month); and • main source of drinking water (8 levels: Piped into dwelling, Piped to yard/plot, Public tap/standpipe, Protected well, Unprotected well, River/dam/lake/ponds/stream/canal/ irrigation channel, Bottled water, Other).

C Discussion: the PDC Class and Other Model Section Criteria
Although we defined the d-class in a general manner for two models F θ1 and F θ2 , most model selection criteria are based on a chosen divergence between a given candidate model F θ and the true model F 0 . In such a setting, we may simplify our definition given in (5) and write instead which leads to the following (consistent) estimator C = D g 1 (y) , g 2 y,θ + tr cov g 1 (y) , ∇ψ g 2 y,θ . (A-2) Mallow's C p can be derived from (A-2). Indeed, consider the linear model given in Setting A. Suppose that we wish to find an estimator for the following criterion: Clearly, this model selection criterion belongs to the d-class error (and to the q-class as well) with g 1 (Y 0 ,β 1 ) = Y 0 , g 2 (Y ,β 2 ) = Xβ and ψ (z) = z T z. By applying Theorem 1 we obtain ∆ = tr cov Y , 2Xβ = 2σ 2 ε tr (S) = 2σ 2 ε p where S denotes the "hat" matrix of X, i.e. S = X X T X −1 X T . Using (A-2) we get C = D g 1 y,θ 1 , g 2 y,θ 2 + ∆ = ||y − Xβ|| 2 2 + 2σ 2 ε p which is unsurprisingly equal to Mallow's C p .
At this point, it must be pointed out that many criteria can be defined using (A-1). However, only a few of them are meaningful for the task of model selection. Indeed, an estimator of C in (A-1) is used in practice to compare several models, and in a stepwise procedure, two models at the time. Hence, C should ideally satisfy (at least) two properties. Consider two candidate models, say F θ1 and F θ2 such that F θ1 is nested within F θ2 , then, C should satisfy the following two properties: (Pr.1) For nested models, the apparent divergence is strictly non-increasing with respect to model complexity. More formally, we have that for any nested models F θ1 in F θ2 the apparent divergence is such that: D g 1 (y) , g 2 y,θ 1 ≥ D g 1 (y) , g 2 y,θ 2 .
(Pr.2) For nested models, the divergence optimism is strictly non-decreasing with respect to model complexity. More formally, we have that for any nested models F θ1 in F θ2 the divergence optimism is such that: tr cov g 1 (y) , ∇ψ g 2 y,θ 1 ≤ tr cov g 1 (y) , ∇ψ g 2 y,θ 2 .
If Property (Pr.1) is not satisfied, it would mean that a larger model (i.e. F θ2 ) can provide a poorer (apparent) fit than F θ1 . Such situations are extremely unlikely. In Corollary A-1 (below) we relate Property (Pr.1) to the estimatorθ. This enables to easily construct a criterion which verifies (Pr.1) given an estimatorθ. The proof of this result is straightforward and therefore shall be omitted.
The rational behind Property (Pr.2) is the following. Suppose that F θ1 is the true model (or is nested within the true model) and assume that models F θj with j = 1, ..., K are models with increasing complexity such that model F θj is nested within model F θ k if 1 ≤ j < k ≤ K. Then, if a criterion satisfies Property (Pr.1), the apparent divergence will be smaller as the model complexity increases. Therefore, the optimism should have exactly the opposite property in order to allow a model "close" to F θ1 to be selected by the criterion at hand.

D.1 Proof of Theorem 1
Proof: Since the divergence D (·, ·) belongs to the d-class we have that D (·, ·) is a valid Bregman divergence, we may express (1) as: So by subtracting the two terms we obtain Let x and z be two real-valued vectors of the same dimension such that E x T x < ∞ and E z T z < ∞. Then we can always write Using (A-3) we can express ∆ as Since E 0 g 1 Y 0 ,θ 1 is a non-stochastic quantity, it follows that which verifies the result of Theorem 1.

D.2 Proof of Theorem 2
To prove Theorem 2, we first state and prove the following Lemma.
Lemma A-1: Suppose that we have three competing models which can be written as: where X (1) , X (2) and X (3) are respectively n×r 1 , n×r 2 and n×r 3 matrices and the dimensions of β 1 , of β 2 and of β 3 are defined accordingly. We assume that the matrix X = X (1) X (2) X (3) is a non-random n × (r 1 + r 2 + r 3 ) full rank matrix. Moreover, under model M l we suppose that ε l ∼ N 0, σ 2 ε I , l = 1, 2, 3. Let alsoŶ l denote the prediction function associated to model M l , l = 1, 2, 3. Then the following results hold under the conditions of Theorem 1: 1. Assuming M 1 or one of its sub-models to be the correct model, we have Assuming M 2 to be the correct model, we have 3. Assuming M 3 to be the correct model, we have Proof: By definition ||Ŷ 1 −Ŷ 2 || 2 2 is a divergence measure. So by using (6) for the squared loss function, we have that and under the conditions of Theorem 1 we have that E Ŷ T jŶ j < ∞, j = 1, 2. So we can write and we obtain E ||Ŷ 1 −Ŷ 2 || 2 2 = tr cov Ŷ 1 ,Ŷ 1 + tr cov Ŷ 2 ,Ŷ 2 − 2 tr cov Ŷ 1 ,Ŷ 2 + tr E Ŷ 1 E T Ŷ 1 Since X is full rank the matrices S (1) and S (2) exist and the first three elements of (A-7) can be simplified under model M l , l = 1, 2, 3, as However, the last three elements of (A-7) depend on the assumed true model. Indeed, if M 1 is the correct model, we have E Ŷ 1 = E Ŷ 2 and, therefore By using (A-8) and (A-9), we have that (A-7) simplifies to and we obtain (A-10)

Using (A-8) and (A-10) simplifies (A-7) to
thus verifying (A-5). Finally, if M 3 is the correct model and since X is full rank, we have (A-11) By using (A-8) and (A-11) we have that (A-7) simplifies to thus verifying (A-6).
Using Lemma A-1, the proof of Theorem 2 is the following.
Proof: Using the definition of the PDC and since σ 2 ε is either known or estimated by an unbiased estimator, we have that it follows, using (A-4), that In order to demonstrate the results of Theorem 2, we have to show that E PDC j,j+m is minimized for j =  0 and m = 1. We consider three different cases: 1. 0 ≤ j + m ≤  0 : In this case the models considered for estimatingβ j andβ j+m can be rewritten as: where X j+1:j+m is a n × m matrix that contains the columns j + 1 to j + m of the matrix X. Under this setting we can use (A-6) of Lemma A-1 and obtain and where X j+m+1:0 is a matrix that contains the columns j + m + 1 to  0 of the matrix X. Compared to (A-13), we have that for sufficiently large n since both B T B and C T C are positive strictly increasing functions of n. Hence for sufficiently large n and j + m ≤  0 , (17) is verified.
2. 0 < j <  0 and  0 < j + m ≤ K + 1: In this case we consider the following models for estimating β j and β j+m : It can be noted that in this case we have that so using (A-7) and (A-15) leads to Moreover, using (A-10) in the proof of Lemma A-1 yields Compared to (A-13), we have that, for sufficiently large n, since B T B is a positive strictly increasing function of n. Hence for sufficiently large n, (17) is also verified when 0 < j <  0 and  0 < j + m ≤ K + 1.

D.3 Proof of Theorem 3
Before considering the proof of Theorem 3, we state and prove the following lemmas.
Using Lemma A-3, the proof of Theorem 3 is the following.

From the proof of Lemma A-3 and by combining (A-23) and (A-24) we get
where Z ∼ N (0, 1) and c is a positive constant. Moreover, by the tail probability of the standard normal distribution we get the following inequality: Thus, we have Then we can write Then we have that b j ≤ a j ≤ c j , ∀j, ∀n since λ n > 0 by (A.2), γ n > 0 by (A-18) in Lemma A-2 and j ≥ 0. Clearly, In addition, we have that Using (A-18) in Lemma A-2 we have that γ 2 n = O (n) and by Assumption (A.2) we have that λ n  0 = O ( √ n) and consequently that  0 = O ( √ n). Thus, we may write We also have that and since the matrices S j+m − S j and S 0+m − S 0 are symmetric and idempotent, we have Craig, 1943, Theorem 2 or Appendix H). Now, let Then, the Moment Generating Function (MGF) associated to X 1 and X 2 is Let Y follow a variance-gamma distribution. From (A-37) we know that the MGF associated to Y is given by Clearly, with the set of parameters µ = 0, β = 0, α = 1 2 and ξ = m 2 the MGFs of Y and G m are the same for all t in some neighborhood of 0. This implies that F Gm (u) = F Y (u) , ∀u. Thus, G m follows a modified variance-gamma distribution whose PDF is given in (A-38). Therefore, when σ 2 ε is known, we have In the case where σ 2 ε is unknown, Assumption (A.1) enables to write So by Slutsky's theorem we obtain that which concludes the proof.

D.5 Asymptotic Probability of Overfitting
Let Z 2 l iid ∼ χ 2 1 , l = 1, ..., k + 1. Then we define the random variable W l as W l = Z 2 l − Z 2 l+1 , l = 1, 2, ..., k . In addition let a ∈ N and b ∈ N be such that k = a + b. We define the following probability: Then, the asymptotic probability of overfitting for PDC λn j,j+1 is given in the following theorem. where 1 denotes the indicator function.
To prove Theorem A-7 we first state and prove the following lemmas.
Proof: Let A and B be two idempotent n×n matrices and let Y ∼ N µ, σ 2 ε I . By Craig (1943, Theorem 2) (see also Appendix H) we have that two quadratic forms, say q 1 = Y T AY and q 2 = Y T BY , are independently distributed as χ 2 if AB = 0. Therefore, all Z 2 l , l = 1, ..., K −  0 + 1 are all mutually independent in the case where S i,i = (S i+1 − S i ) (S i +1 − S i ) = 0 for any i ∈ J 0 , i ∈ J 0 and i > i . Consider the term S i,i . We have that and thus we verify that all Z 2 l , l = 1, ..., K −  0 + 1 are all mutually independent and conclude the proof.
The proof of Theorem A-7 can now be given.
Proof: This proof is straightforward from Lemma A-5 and Theorem 3. Nevertheless, we provide here a proof for completeness. By applying the law of total probability we have the following equality: Since by Assumption (A.1) we have that then by Slutsky's theorem, which concludes the proof.

D.6 Proof of Theorem 5
To prove Theorem 5, we first need to state and prove the corrolaries.
Corollary A-2: Let µ = Xβ 0 + δ. Then if J = J 0 and if σ 2 ε is known, we have that and under the assumptions of Theorem A-7, we have that Since J = J 0 we have that From Theorem A-7 we have the following Corollary Corollary A-3: Under the assumptions of Theorem A-7 we have that Proof: The proof is straightforward from Theorem A-7: Remark G: In the situation where λ n = 2, we obtain from Corollary A-3 that while for the C p (and the AIC) we have that (see Woodroofe, 1982 andZhang, 1992 for details).
From Corollary A-2 and Corollary A-3, we may define the asymptotic efficiency as In the case where δ = 0, we simply obtain Clearly if λ n is such that it tends to infinity with n, the asymptotic efficiency tends to 1.
Remark H: In the case where λ n = 2, we have approximately that which is approximately equal to 1 for large values of  0 .
Proof of Theorem 5: Note that some elements of this theorem are closely related to Theorem 5.3 of Müller (2008).
First, we define for notational simplicity B n (j) as Second, by applying the law of total probability, we have for any given > 0 that The result of Theorem 3 implies that lim n→∞ Pr ( λn / ∈ J 0 ) = 0. So we only need to show that lim n→∞ Pr (|B n ( λn ) − 1| ≥ | λn ∈ J 0 ) = 0 to prove that B n ( λn ) P − → 1: where a j = L n ( 0 ) − υ n and where υ = µ T S λn − S 0 µ. Note that using (A-19) in Lemma A-2 we have that υ = O (n). Using Markov's inequality yields to Finally, by taking the limit we can write where c is a positive constant. Therefore, we have that B n ( λn ) P − → 1.
Next, we consider the case where lim n→∞ λ n = ∞. In such situations, we have that

D.7 Proof of Theorem 6
Proof: At step j, the variable indexed by k will be selected if Let K j = {1, ..., K − j} and K j denotes the set of indices in K j corresponding to the K − j unordered variables such that these variables are all associated to one of the first  0 columns of the matrix X . In order to verify Theorem 6 we have to prove that at step j, and assuming that K j is not empty, k asymptotically belongs to K j .
Clearly, we have that j+1 denotes the "hat" matrix associated toŶ (k ) j+1 . Note that by (A-18) in Lemma A-2 we have that γ 2 n = O (n). Let Y 1,γn ∼ χ 2 1 γ 2 n and Z 1 ∼ N (0, 1). We can always write that Y 1,γn = (Z 1 + γ n ) j+1 || 2 2 also follows a χ 2 1 distribution and is independent from σ −2 ||Ŷ j −Ŷ (k ) j+1 || 2 2 by Craig (1943, Theorem 2) (see also Appendix H) and therefore is also independent from Z 1 . Let Z 2 ∼ N (0, 1) be independent from Z 1 and let By means of Chebyshev's inequality we can show that U = o p (1) (see Lemma A-3). Moreover, we also have that where c is a positive constant since γ n > 0 and γ n = O ( √ n). By Slutsky's theorem we have We can rewrite W as follows: Finally, since n − 1 /2 W = o p (1) and we also have that where c is a positive constant. Therefore, we have that, asymptotically, the probability that k belongs to K j is one, which proves Theorem 6.

E Signal to noise ratio Theorem
Suppose that one wishes to compare two nested models M 1 (with k variables) and M 2 (with m ≥ 1 additional variables) using an estimator of the criterion C. We define the ∆ C as where C Mj , j = 1, 2 denotes the values of C for model M j . Clearly, model M 1 will be considered better in the case where ∆ C ≤ 0 while if ∆ C > 0 model M 2 will be considered preferable. The random variable ∆ C depends largely on the strength of the penalty function (i.e. the divergence optimism) and the SNR is a possible tool that can be used to study the decision rule (A-27). Following McQuarrie and Tsai (1998) we define the SNR as where the numerator of the SNR corresponds to the signal and the denominator to the noise. As explained in McQuarrie and Tsai (1998), while the signal depends primarily on the difference of divergence optimism of the criterion used, the noise depends mainly on the distribution of the difference of apparent errors. A model selection procedure with a small penalty function will tend to have a weak signal, a weak SNR and will be prone to overfit. Although a procedure with a larger SNR might overcome this difficulty, such an approach might tend to underfit. The AIC is an example of the former while the BIC is one of the latter. In the following theorem, the SNR of PDC is compared to the SNR of some of the model selection criteria presented in Table A-5 in Appendix F, assuming σ 2 ε known. More general results (in the overfitting case) can be found in McQuarrie and Tsai (1998, Section 2.3).
The AIC, BIC and HQ are based on log-likelhood, i.e n log (2π) + n log σ 2 ε + 1 and can be expressed as Crit = ||Y −Ŷ || 2 2 + σ 2 ε penalty + C where C does not depend on the candidate model and can therefore be neglected.
Since the penalty term of the AIC is 2p, it is equivalent to the C p and therefore it has the same SNR, thus verifying (A-30).
For the BIC, where the penalty term is log(n)p. We can write By using (A-34) we obtain, Therefore, the SNR of the BIC is given by which verifies (A-31).
Therefore, the SNR of the HQ is given by which verifies (A-32).
Finally, we obtain which verifies (A-33) and concludes the proof.
The left panel of Figure  as a function of β. Globally, it can be observed that the PDC has a larger SNR than the C p and the same can be said for PDC and PDC • with respect to the HQ and the BIC. As it can be seen in the right panel of Figure A-5 that a criterion with a high SNR tends to choose the smallest model M 1 more frequently. Therefore, criteria of the PDC class select on average smaller models than the "classical" counterpart. In this example, the PDC has a behavior that is very similar to the HQ information criterion and is a compromise between the behavior of the BIC and C p .

F Criteria for model selection in linear regression
In Table A-5 are provided some of the most commonly used criteria for model selection in the context of linear regression models. For the lasso and adaptive methods, we used the following specifications: • For the lasso we used the R function lars of the lars package. The shrinkage coefficient λ was chosen by minimizing the C p statistic.
• For the elastic-net we used the R function enet of the elasticnet package. The shrinkage coefficient λ 2 (L 2 penalty) was chosen by tenfold cross validation and λ 1 (L 1 penalty) by minimizing the C p statistic.
• For the adaptive lasso we used the R function adalasso of the parcor which computes the solution based on tenfold cross validation. The initial weights for adaptive lasso are computed from a lasso fit. In the special case where m = 1, (A-38) can be simplified as Note that the function K 0 (z) has the following integral representation (see e.g. Abramowitz and Stegun, 1965, eq. 9.6.22, p. 376):

H Basic results on Quadratic Forms
Let Y ∼ N µ, σ 2 I be an n-vector with a spherical normal distribution and A be an n × n symmetric idempotent matrix with rank (A) = r. Then A proof of (A-41) can, for example, be found in Bates (2011, proof

I Simulation results
We present here the results of the simulation studies 5.1, 5.2 and 5.3. The evaluation criteria are presented in Table 3 and the simulation results are reported in Tables A-6 to A-8.  Table 3. These include the PDC λn j,j+1 criteria in (16), with λn = 2 (P DC), λn = log(n) (PDC • ), λn = log(log(n)) (PDC ), the stepwise forward FPE (FPE), FPEu (FPEu), AIC (AIC), AICc (AICc), AICu (AICu), BIC (BIC), HQ (HQ), HQc (HQc), and the full model with the LSE ( LSE), together with the lasso (lasso), elastic-net (enet), adaptive lasso (a-lasso), MCP (M CP ) and SCAD (SCAD) (see Appendix F for more details on the chosen criteria). 500 samples are simulated under the correct model as presented in Simulation 5.1 (p = 2, K = 60). The numbers in parentheses for the columns Med (PEy) and Med (MSE β ) are the corresponding standard errors estimated by using the bootstrap with B = 500 resamplings. The numbers in superscript indicate the (average) ranked performance for each evaluation criterion (before rounding).