Selection consistency of Lasso-based procedures for misspecified high-dimensional binary model and random regressors

TL;DR

This paper studies the consistency of Lasso-based methods for selecting predictors in high-dimensional binary regression models, including misspecified and semiparametric cases, demonstrating their effectiveness under certain conditions.

Contribution

It introduces a two-step Lasso-based selection procedure with proven consistency in high-dimensional, potentially misspecified binary models, extending existing theoretical guarantees.

Findings

The proposed method is consistent for large predictor sets.

Support of true parameters can be accurately recovered.

Method applies to semiparametric models with linear predictor conditions.

Abstract

We consider selection of random predictors for high-dimensional regression problem with binary response for a general loss function. Important special case is when the binary model is semiparametric and the response function is misspecified under parametric model fit. Selection for such a scenario aims at recovering the support of the minimizer of the associated risk with large probability. We propose a two-step selection procedure which consists of screening and ordering predictors by Lasso method and then selecting a subset of predictors which minimizes Generalized Information Criterion on the corresponding nested family of models. We prove consistency of the selection method under conditions which allow for much larger number of predictors than number of observations. For the semiparametric case when distribution of random predictors satisfies linear regression conditions the true and the estimated parameters are collinear and their common support can be consistently identified.

Equations

$$E\exp (t{X}_{ij})\le \exp (t^2{\sigma }_{jn}^2/2)$$

$$s_n^2=\underset{j=1,\dots ,p_n}{\max }{\sigma }_{jn}^2$$

$${\gamma }^2:=\underset{n}{\mathrm{limsup}}{s}_n^2<\infty .$$

$$l(b,x,y)=\rho (b^{T}x,y),$$

$$R(b)=El(b,x,y)$$

$${\beta }^{\ast }={\mathrm{argmin}}_{b\in R^{p_n}}R(b).$$

$$q(x)=:q({\beta }^{T}x)$$

$$s^{\ast }=\mathrm{supp}({\beta }^{\ast }(\{1,\dots ,p_n\})=\{i\le p_n:{\beta }_i^{\ast }\ne 0\},$$

$${\beta }_{s^{\ast }}^{\ast }={\beta }^{\ast }(s^{\ast }{)}_{s^{\ast }}={\beta }^{\ast }(s{)}_{s^{\ast }}.$$

$${\beta }_{min}^{\ast }=\underset{i\in s^{\ast }}{\min }|{\beta }_i^{\ast }|.$$

$${\hat{\beta }}_L={\mathrm{argmin}}_{b\in R^{p_n}}\{R_n(b)+{\lambda }_L\sum _{i=1}^{p_n}|b_i|\}$$

$$R_n(b)=\sum _{i=1}^n\rho (b^T{X}_i,Y_i).$$

$${\mathcal{C}}_{\epsilon }=\{\Delta \in R^{p_n}:\text{leavevmode}||{\Delta }_{s^{\ast c}}|{|}_1\le (3+\epsilon )||{\Delta }_{s^{\ast }}|{|}_1\},$$

$${\kappa }_H(\epsilon )=\underset{\Delta \in {\mathcal{C}}_{\epsilon }\setminus \{0\}}{\inf }\frac{{\Delta }^{T}H\Delta }{{\Delta }^T\Delta },$$

$$W(b)=R(b)-R({\beta }^{\ast }),$$

$$W_n(b)=R_n(b)-R_n({\beta }^{\ast }),$$

$$B_p(r)=\{\Delta \in R^{p_n}:\text{leavevmode}||\Delta |{|}_p\le r\},\text{ for }p=1,2,$$

$${\beta }_{min}^{\ast }=\underset{i\in s^{\ast }}{\min }|{\beta }_i^{\ast }|,$$

$$S(r)=\underset{b\in R^{p_n}:b-{\beta }^{\ast }\in B_1(r)}{\sup }|W(b)-W_n(b)|.$$

$$R(b)-R({\beta }^{\ast })\ge \frac{\vartheta }{2}(b-{\beta }^{\ast }{)}^{T}H(b-{\beta }^{\ast }).$$

$$R_n({\beta }^{\ast }+\Delta )-R_n({\beta }^{\ast })\ge D{R}_n({\beta }^{\ast }{)}^T\Delta +{\kappa }_L||\Delta |{|}^2-{\tau }^2({\beta }^{\ast })$$

$$GIC(w)=n{R}_n(\hat{\beta }(w))+a_n|w|,$$

$$\hat{\beta }(w)=\underset{b\in R^{p_n}:\text{leavevmode}{b}_{w^c}=0_{|w^c|}}{\mathrm{argmin}}{R}_n(b)$$

$$S_k(r)=\underset{b\in D_k:b-{\beta }^{\ast }\in B_2(r)}{\sup }|(W_n(b)-W(b)|$$

$$D_1=\{b\in R^{p_n}:\text{leavevmode}\exists w\in \mathcal{M}:\text{leavevmode}|w|\le k_n\wedge s^{\ast }\subset w\wedge \mathrm{supp}b\subseteq w\},$$

$$D_2=\{b\in R^{p_n}:\text{leavevmode}\mathrm{supp}b\subset s^{\ast }\}.$$

$$R(b)-R({\beta }^{\ast })\ge \frac{\vartheta }{2}(b-{\beta }^{\ast }{)}^{T}H(b-{\beta }^{\ast })\ge \frac{\vartheta {\lambda }_{min}}{2}||b-{\beta }^{\ast }|{|}_2^2.$$

$$R(b)-R({\beta }^{\ast })\ge \theta ||b-{\beta }^{\ast }|{|}_2^2\ge \frac{\theta }{{\lambda }_{max}}(b-{\beta }^{\ast }{)}^{T}H(b-{\beta }^{\ast }).$$

$$u=\frac{r}{r+||{\hat{\beta }}_L-\beta |{|}_1},v=u{\hat{\beta }}_L+(1-u){\beta }^{\ast },$$

$$W(v)+\lambda ||v-{\beta }^{\ast }|{|}_1\le S(r)+2\lambda ||v_{s^{\ast }}-{\beta }_{s^{\ast }}^{\ast }|{|}_1.$$

$$P(S(r)>t)$$

$$\le \frac{2}{tn}E\underset{b\in R^{p_n}:b-{\beta }^{\ast }\in B_1(r)}{\sup }\left|\sum _{i=1}^n{\epsilon }_i(\rho (X_i^{T}b,Y_i)-\rho (X_i^T{\beta }^{\ast },Y_i))\right|$$

$$\le \frac{4L}{tn}E\underset{b\in R^{p_n}:b-{\beta }^{\ast }\in B_1(r)}{\sup }\left|\sum _{i=1}^n{\epsilon }_i{X}_i^T(b-{\beta }^{\ast })\right|.$$

$$E{\left|\left|\sum _{i=1}^n{\epsilon }_i{X}_{ij}\right|\right|}_{\infty }\le \sqrt{n}{s}_n\sqrt{2\ln (2p_n)}\le 2s_n\sqrt{n\ln (p_n\vee 2)}$$

Full text

**Keywords: ** high-dimensional regression, loss function, random predictors, misspecification, consistent selection, subgaussianity, Generalized Information Criterion

1 Introduction

We consider random variable $(X^{(n)},Y^{(n)})\in R^{p_{n}}\times\{0,1\}$ and corresponding response function defined as a posteriori probability $q_{n}(x)=P(Y^{(n)}=1|X^{(n)}=x)$. We adopt triangular scenario and assume that $n$ copies $X_{1}^{(n)},\ldots,X_{n}^{(n)}$ of a random vector $X^{(n)}$ in $R^{p_{n}}$ are observed together with corresponding binary responses $Y_{1}^{(n)},\ldots,Y_{n}^{(n)}$. We assume that observations $(X_{i}^{(n)},Y_{i}^{(n)}),\,i=1,\ldots,n$ are iid. Let $X_{i}^{(n)}=(X_{i1}^{(n)},\ldots,X_{ip_{n}}^{(n)})^{\prime}$. Frequently considered scenario is the sequential one. In this case, when sample size $n$ increases we observe values of new predictors additionally to the ones observed earlier. This is a special case of the above scheme as then $X_{i}^{(n+1)}=({X_{i}^{(n)T}},X_{i,p_{n}+1},\ldots,X_{i,p_{n+1}})^{T}$. In the following we will skip the upper index $n$ if no ambiguity arises. Moreover, we write $q(x)=q_{n}(x)$. We assume that coordinates $X_{ij}$ of $X_{i}$ are subgaussian $Subg(\sigma_{jn}^{2})$ with subgaussianity parameter $\sigma_{jn}^{2}$ i.e. it holds that

[TABLE]

for all $t\in R$. For future reference let

[TABLE]

and assume in the following that

[TABLE]

In the sequential scenario this is equivalent to an assumption that all subgaussianity parameters are bounded from above. We assume moreover that $X_{i1},\ldots,X_{i{p_{n}}}$ are linearly independent in the sense that their arbitrary linear combination is not constant almost everywhere. In the following ${\mathbb{X}}_{n}=(X_{1},\ldots,X_{n})^{T}$ will denote matrix of experiment of dimension $n\times p_{n}$.

For the regression defined above we consider loss function of the form

[TABLE]

where $\rho:R\times\{0,1\}\to R$ is some function, $b,x\in R^{p_{n}},$ $y\in\{0,1\}$ and

[TABLE]

is associated risk function for $b\in R^{p_{n}}$. Our aim is to determine the support of $\beta^{*}$, where

[TABLE]

Coordinates of $\beta^{*}$ corresponding to non-zero coefficients will be called active predictors and vector $\beta^{*}$ a pseudo-true vector. This terminology stems from the important special case of our general setting: misspecification case of semiparametric model. Namely, consider a semiparametric model for which response function is given in semiparametric form

[TABLE]

for some fixed $\beta$ and unknown $q$. When the loss defined in (2) does not coincide with minus conditional log-likelihood -$E(Y\log q(b^{T}{X})+(1-Y)\log(1-q(b^{T}{X}))$ pertaining to $q(b^{T}{x})$, in particular when fitted parametric model is given by a response function $q_{0}\not\equiv q$ then the model is misspecified. Questions of robustness analysis evolve around the interplay between $\beta$ and $\beta^{*}$, in particular under what conditions the directions of $\beta$ and $\beta^{*}$ coincide (cf important contribution in Brillinger (1982) and Ruud (1983)).

In the paper we consider properties of $\beta^{*}$ for a general loss function and the case of misspecified semiparametric model (4) as a special case of this setup. For $s\subseteq\{1,\ldots,p_{n}\}$ let $\beta^{*}(s)$ be defined as in (3) when minimum is taken over $b$ with support in $s$. We define

[TABLE]

denote the support of $\beta^{*}(\{1,\ldots,p_{n}\})$ with $\beta^{*}(\{1,\ldots,p_{n}\})=(\beta_{1}^{*},\ldots,\beta_{p_{n}}^{*})^{T}$.

Let $v_{\pi}=(v_{j_{1}},\ldots,v_{j_{k}})^{T}\in R^{|\pi|}$ for $v\in R^{p_{n}}$ and $\pi=\{j_{1},\ldots,j_{k}\}\subseteq\{1,\ldots,p_{n}\}$. Let $\beta^{*}_{s^{*}}\in R^{|s^{*}|}$ be $\beta^{*}=\beta^{*}(\{1,\ldots,p_{n}\})$ restricted to its support $s^{*}$. Note that if $s^{*}\subseteq s$, then provided projections are unique (see Section 2) we have

[TABLE]

Moreover, let

[TABLE]

We remark that $\beta^{*}$, $s^{*}$ and $\beta_{min}^{*}$ may depend on $n$. Note that when the parametric model is correctly specified i.e. $q(x)=q(\beta^{T}x)$ for some $\beta$ with $l$ being an associated loglikelihhood loss and if $s$ is the support of $\beta$ then we have $s=s^{*}$.

For fixed number $p$ of predictors smaller than sample size $n$ statistical consequences of misspecification of a semiparametric regression model have been intensively studied by H. White and his collaborators in the 80s of the last century. The concept of projection on the fitted parametric model is central to this investigations which show how the distribution of maximum likelihood estimator of $\beta^{*}$ centered by $\beta^{*}$ changes under misspecification (cf e.g. White (1982) and Vuong (1989)). However for the case when $n>p$ the maximum likelihood estimator which is a natural tool for fixed $p<n$ case is ill-defined and a natural question arises what can be estimated and by what means in this case.

For high-dimensional case one of possible solutions is to consider two-stage methods in which the first stage results in screened subset of regressors with cardinality smaller than $n$ and the second stage employs one of known methods for fixed $p$ case. As the set of regressors for the second stage is random the properties of the procedure need to be thoroughly reconsidered.

In the paper first stage of the procedure is is based on Lasso estimation

[TABLE]

where $b=(b_{1},\ldots,b_{p_{n}})^{T}$ and the empirical risk $R_{n}(b)$ is

[TABLE]

Parameter $\lambda_{L}>0$ is Lasso penalty which penalizes large $l_{1}$-norms of potential candidates for a solution. Note that criterion function in (5) for $\rho(s,y)=\log(1+\exp(-s(2y-1))$ can be viewed as penalized empirical risk for the logistic loss. Lasso estimator is thoroughly studied in the case of the linear model when considered loss is square loss (see e.g. Bühlmann and van de Geer (2011) and Hastie et al. (2015) for references and overview of the subject) and some of the papers treat the case when such model is fitted to $Y$ which is not necessarily linearly dependent on regressors (cf Bickel et al. (2009) ). In this case regression model is misspecified w.r.t. linear fit. However, similar results are scarse for other scenarios such as logistic fit under misspecification in particular. One of the notable exceptions is Negahban et al. (2012) where behaviour of Lasso estimate is studied for a general loss function including logistic loss for possibly misspecified models. For a recent contributions to study of Kullback-Leibler projections on logistic model (which coincide with (3) for logistic loss) and references we refer to Kubkowski and Mielniczuk (2017) and Kubkowski and Mielniczuk (2018). We also refer to Lu et al. (2012), where asymptotic distribution of adaptive Lasso is studied under misspecification in the case of fixed number of deterministic predictors. In the following we prove approximation result for Lasso estimator when predictors are random and $\rho$ is a convex Lipschitz function (cf Theorem 1). An useful corollary of it is determination of sufficient conditions under which active predictors can be separated from spurious ones based on the absolute values of corresponding coordinates of Lasso estimator. This makes construction of nested family containing $s^{*}$ with large probability possible. In the general framework allowing for misspecification we call selection rule $\hat{s}^{*}$ consistent if $P(\hat{s}^{*}=s^{*})\to 1$ when $n\to\infty$.

In the case of semiparametric model when predictors are eliptically contoured ( e.g. multivariate normal) it is known that $\beta^{*}=\eta\beta$ i.e. these two vectors are collinear (Li and Duan (1989)). Thus in case when $\eta\neq 0$ we have that $s^{*}$ coincides with support $s$ of $\beta$ and the selection consistency of two-step procedure proved in the paper entails direction and support recovery of $\beta$.

The main objective of the paper is to prove consistency of two-stage selection procedure which consists of ordering of predictors according to the absolute values of corresponding Lasso estimators and then minimization of Generalized Information Criterion GIC) on resulting nested family.This is a variant of SOS (Screening-Ordering-Selection) procedure introduced in Pokarowski and Mielniczuk (2015) in the case of the linear model, where the ordering of predictors chosen by Lasso was with respect to absolute values of $t$ statistics from the linear fit based on these predictors. Here we consider a simpler scheme for which both screening and ordering of regressors is based on Lasso fit. We consider the case when predictors are subgaussian random variables. The stated results to the best of our knowledge are not available for random predictors even when the model is correctly specified. For the second stage we assume that the number of active predictors is bounded by a deterministic sequence $k_{n}$ tending to infinity and we minimize GIC on family ${\cal M}$ of models with sizes satisfying also this condition. Such exhaustive search has been proposed in Chen and Chen (2008) for linear models and extended to GLMs in Chen and Chen (2012), see also Mielniczuk and Szymanowski (2015). In these papers GIC has been optimised on all possible subsets of regressors with cardinality not exceeding certain constant $k_{n}$. Such method is feasible for practical purposes only when $p_{n}$ is small. Here we consider a similar setup but with important differences: ${\cal M}$ is a data-dependent small nested family of models and optimization of GIC is considered in the case when the original model is misspecified. The regressors are assumed random and assumptions are carefully tailored to this case.

In numerical experiments we study the performance of grid version of logistic and linear SOS and compare it to its several lasso-based competitors.

The paper is organized as follows. Section 2 contains auxiliaries, including new useful probability inequalities for empirical risk in the case of subgaussian random variables (Lemma 2). In Section 3 we prove a bound on approximation error for Lasso for misspecified logistic model and random regressors (Theorem 1) which yields separation property of Lasso. In Theorems 2 and 3 of Section 4 we prove GIC consistency on nested family, which in particular can be built according to the order of the Lasso coordinates. In Corollary 5 we discuss consequences of the proved semiparametric binary model when distribution of predictors satisfies linear regressions condition. In Section 5 we compare the performance of two-stage selection method for two closely related models, one of which is a logistic model and the second one is misspecified.

2 Definitions and auxiliary results

We assume throughout existence and uniqueness of projection vector $\beta^{*}$ which has been defined in (3). We consider cones of the form:

[TABLE]

where $\varepsilon>0$, $s^{*c}=\{1,\ldots,p_{n}\}\setminus s^{*}$ and $\Delta_{s^{*}}=(\Delta_{s^{*}_{1}},\ldots,\Delta_{s^{*}_{|s^{*}|}})$ for $s^{*}=\{s^{*}_{1},\ldots,s^{*}_{|s^{*}|}\}$. Cones $\mathcal{C}_{\varepsilon}$ are of special importance (see Lemma 3). For cone $\mathcal{C}_{\varepsilon}$ we define a quantity $\kappa_{H}(\varepsilon)$ which can be regarded as generalized minimal eigenvalue of a matrix in high-dimensional setup:

[TABLE]

where $H\in R^{p_{n}\times p_{n}}$ is non-negative definite matrix, which in the considered context is usually taken as hessian $D^{2}R(\beta^{*})$ .

Let $R$ and $R_{n}$ be the risk and the empirical risk defined above. Moreover, we introduce the following notation:

[TABLE]

We will need the following margin condition in Lemma 3 and Theorem 1:

1. (MC)

There exist $\vartheta,\varepsilon,\delta>0$ and non-negative definite matrix $H\in R^{p_{n}\times p_{n}}$ such that for all $b$ with $b-\beta^{*}\in\mathcal{C}_{\varepsilon}\cap B_{1}(\delta)$ we have

[TABLE]

The above condition can be viewed as a weaker version of strong convexity of function $R$ in the restricted neighbourhood of $\beta^{*}$ (namely in the intersection of ball $B_{1}(\delta)$ and cone $\mathcal{C}_{\varepsilon}$). We stress the fact that $H$ does not need to be positive definite, as in the Section 3 we use (MC) together with stronger conditions than $\kappa_{H}(\varepsilon)>0$ which imply that right hand side of inequality in (MC) is positive. We also do not require here twice differentiability of $R$. We note in particular that condition (MC) is satisfied in the case of logistic loss, $X$ being bounded random variable and $H=D^{2}R(\beta^{*})$ (see Fan et al. (2014b) and Bach (2010)). It is also easily seen that that (MC) is satisfied for quadratic loss, $X$ satisfying $E||X||_{2}^{2}<\infty$ and $H=D^{2}R(\beta^{*})$. Similar condition to (MC) (called Restricted Strict Convexity) was considered in Negahban et al. (2012) for empirical risk $R_{n}$:

[TABLE]

for all $\Delta\in C(3,s^{*})$, some $\kappa_{L}>0$ and tolerance function $\tau$.

Another important assumption, used in the Theorem 1 and Lemma 2 is the Lipschitz property of $\rho:$

1. (LL)

$\exists L>0\ \forall b_{1},b_{2}\in R,y\in\{0,1\}\colon\ |\rho(b_{1},y)-\rho(b_{2},y)|\leq L|b_{1}-b_{2}|$.

Let $|w|$ stand for dimension of $w$. For the second step of the procedure we consider an arbitrary family $\mathcal{M}\subseteq 2^{\{1,\ldots,p_{n}\}}$ of models (which are identified with subsets of $\{1,\ldots,p_{n}\}$ and may be data-dependent) such that $s^{*}\in\mathcal{M},\forall w\in\mathcal{M}:\ |w|\leq k_{n}$ a.e. and $k_{n}\in N_{+}$ is some deterministic sequence. We define Generalized Information Criterion (GIC) as:

[TABLE]

where

[TABLE]

is ML estimator for model $w$ and $a_{n}>0$ is some penalty. Typical examples of $a_{n}$ include:

• •

AIC (Akaike Information Criterion): $a_{n}=2$,

• •

BIC (Bayesian Information Criterion): $a_{n}=\log n$,

• •

EBIC($d$) (Extended BIC): $a_{n}=\log n+2d\log p_{n}$, where $d>0$.

We will study properties of $S_{k}(r)$ for $k=1,2$, where:

[TABLE]

and is the maximal absolute value of the centred empirical risk $W_{n}(\cdot)$ and sets $D_{k}$ for $k=1,2$ are defined as follows:

[TABLE]

We note that such definitions of $D_{i}$ for $i=1,2$ guarantee that if $b\in D_{i}$, then $|\operatorname*{supp}(b-\beta^{*})|\leq k_{n}$, what we exploit in Lemma 2. Moreover, in Section 4 we consider the following condition for $\epsilon>0$, $w\subseteq\{1,\ldots,p_{n}\}$ and some $\theta>0$:

1. $C_{\epsilon}(w)\colon$

$R(b)-R(\beta^{*})\geq\theta||b-\beta^{*}||_{2}^{2}$ for all $b\in R^{p_{n}}$ such that $\operatorname*{supp}b\subseteq w$ and $b-\beta^{*}\in B_{2}(\epsilon).$

We observe also that the conditions (MC) and $C_{\epsilon}(w)$ are not equivalent, as they hold for $v=b-\beta^{*}$ belonging to different sets: $B_{1}(r)\cap\mathcal{C}_{\varepsilon}$ and $B_{2}(\epsilon)\cap\{\Delta\in R^{p_{n}}\colon\leavevmode\nobreak\ \operatorname*{supp}\Delta\subseteq w\}$, respectively. We note that if the minimal eigenvalue $\lambda_{min}$ of matrix $H$ in condition (MC) is positive and (MC) holds for $b-\beta^{*}\in B_{1}(r)$ (instead of for $b-\beta^{*}\in\mathcal{C}_{\varepsilon}\cap B_{1}(r)$) then we have for $b-\beta^{*}\in B_{2}(r/\sqrt{p_{n}})\subseteq B_{1}(r)$:

[TABLE]

Furthermore, if $\lambda_{max}$ is the maximal eigenvalue of $H$ and $C_{\epsilon}(w)$ holds for all $v=b-\beta^{*}\in B_{2}(r)$ without restriction on $\operatorname*{supp}b$, then we have for $b-\beta^{*}\in B_{1}(r)\subseteq B_{2}(r)$:

[TABLE]

Similar condition to $C_{\epsilon}(w)$ for empirical risk $R_{n}$ was considered in (Kim and Jeon, 2016, formula (2.1)) in the context of GIC minimization.

It turns out that condition $C_{\epsilon}(w)$ together with $\rho(\cdot,y)$ being convex for all $y$ and satisfying Lipschitz condition (LL) are sufficient to establish bounds which ensure GIC consistency for $k_{n}\ln p_{n}=o(n)$ and $k_{n}\ln p_{n}=o(a_{n})$ (see Corollaries 2 and 3).

Lemma 1**.**

(Basic inequality). Let $\rho(\cdot,y)$ be convex function for all $y.$ If for some $r>0$ we have:

[TABLE]

then:

[TABLE]

Proof of the lemma is moved to Appendix. Quantities $S_{k}(r)$ are diefined in (14).

Lemma 2**.**

Let $\rho(\cdot,y)$ be convex function for all $y$ and satisfy Lipschitz condition (LL). Assume that $X_{ij}$ for $j\geq 1$ are subgaussian $Subg(\sigma_{jn}^{2})$, where $\sigma_{jn}\leq s_{n}$. Then for $r,t>0$:

$P(S(r)>t)\leq\frac{8Lrs_{n}\sqrt{\log(p_{n}\vee 2)}}{t\sqrt{n}}$, 2. 2.

$P(S_{1}(r)\geq t)\leq\frac{8Lrs_{n}\sqrt{k_{n}\ln(p_{n}\vee 2)}}{t\sqrt{n}}$, 3. 3.

$P(S_{2}(r)\geq t)\leq\frac{4Lrs_{n}\sqrt{|s^{*}|}}{t\sqrt{n}}$.

Proof.

From the Chebyshev inequality (first inequality below), symmetrization inequality (see (van der Vaart and Wellner, 1996, Lemma 2.3.1)) and Talagrand - Ledoux inequality ((Ledoux and Talagrand, 1991, Theorem 4.12)) we have for $t>0$ and $(\varepsilon_{i})_{i=1,\ldots,n}$ being Rademacher variables independent of $(X_{i})_{i=1,\ldots,n}$:

[TABLE]

We observe that $\varepsilon_{i}X_{ij}\sim Subg(\sigma_{jn}^{2}).$ Hence, using independence we obtain $\sum\limits_{i=1}^{n}\varepsilon_{i}X_{ij}\sim Subg(n\sigma_{jn}^{2})$ and thus $\sum\limits_{i=1}^{n}\varepsilon_{i}X_{ij}\sim Subg(ns_{n}^{2}).$ Applying Hölder inequality and the following inequality (see (Devroye and Lugosi, 2012, Lemma 2.2)):

[TABLE]

we have:

[TABLE]

From this part 1 follows. In the proofs of parts 2-3 first inequalities are identical as in (17) with supremums taken on corresponding sets. Using Cauchy-Schwarz inequality, inequality $||v||_{2}\leq\sqrt{||v||_{0}}||v||_{\infty}$, inequality $||v_{\pi}||_{\infty}\leq||v||_{\infty}$ for $\pi\subseteq\{1,\ldots,p_{n}\}$ and (18) yields:

[TABLE]

Similarly for $S_{2}(r)$, using Cauchy-Schwarz inequality, $||v_{\pi}||_{2}\leq||v_{s^{*}}||_{2}$ which is valid for $\pi\subseteq s^{*}$, definition of $l_{2}$ norm and inequality $E|Z|\leq\sqrt{EZ^{2}}\leq\sigma$ for $Z\sim Subg(\sigma^{2})$, we obtain:

[TABLE]

∎

3 Properties of Lasso for a general loss function and random predictors

The main Theorem in this section is Theorem 1. Idea of the proof is based on fact that if $S(r)$ defined in (12) is sufficiently small, then $\hat{\beta}_{L}$ lies in a ball $\{\Delta\in R^{p_{n}}\colon\leavevmode\nobreak\ ||\Delta-\beta^{*}||_{1}\leq r\}$ (see Lemma 3). Using a tail inequality for $S(r)$ proved in Lemma 2 we obtain Theorem 1. Convexity of $\rho(\cdot,y)$ below is understood as convexity for both $y=0,1$.

Lemma 3**.**

Let $\rho(\cdot,y)$ be convex function and assume that $\lambda>0.$ Moreover, assume margin condition (MC) with constants $\vartheta,\epsilon,\delta>0$ and some non-negative definite matrix $H\in R^{p_{n}\times p_{n}}$. If for some $r\in(0,\delta]$ we have $S(r)\leq\bar{C}\lambda r$ and $2|s^{*}|\lambda\leq\kappa_{H}(\varepsilon)\vartheta\tilde{C}r$, where $\bar{C}=\varepsilon/(8+2\varepsilon)$ and $\tilde{C}=2/(4+\varepsilon),$ then

[TABLE]

Proof.

Let $u$ and $v$ be defined as in Lemma 1. Observe that $||v-\beta^{*}||_{1}\leq r/2$ is equivalent to $||\hat{\beta}_{L}-\beta^{*}||_{1}\leq r,$ as the function $f(x)=rx/(x+r)$ is increasing, $f(r)=r/2$ and $f(||\hat{\beta}_{L}-\beta^{*}||_{1})=||v-\beta^{*}||_{1}$. Let $C=1/(4+\varepsilon).$ We consider two cases:

(i) $||v_{s^{*}}-\beta^{*}_{s^{*}}||_{1}\leq Cr$:

In this case from the basic inequality (Lemma 1) we have:

[TABLE]

(ii) $||v_{s^{*}}-\beta^{*}_{s^{*}}||_{1}>Cr$:

Note that $||v_{s^{*c}}||_{1}<(1-C)r,$ otherwise we would have $||v-\beta^{*}||_{1}>r$ which contradicts (33) in proof of Lemma 1 (see Appendix).

Now we observe that $v-\beta^{*}\in\mathcal{C}_{\varepsilon},$ as we have from definition of $C$ and assumption for this case:

[TABLE]

By inequality between $l_{1}$ and $l_{2}$ norms, definition of $\kappa_{H}(\varepsilon),$ inequality $ca^{2}/4+b^{2}/c\geq ab$ and margin condition (MC) (which holds because $v-\beta^{*}\in B_{1}(r)\subseteq B_{1}(\delta)$ in view of (33)) we conclude that:

[TABLE]

Hence from the basic inequality (Lemma 1) and inequality above it follows that:

[TABLE]

Subtracting $W(v)$ from both sides of above inequality, using assumption on $S,$ the bound on $|s^{*}|$ and definition of $\tilde{C}$ yields:

[TABLE]

∎

Theorem 1**.**

Let $\rho(\cdot,y)$ be convex function for all $y$ and satisfy Lipschitz condition (LL). Assume that $X_{ij}\sim Subg(\sigma_{jn}^{2})$, $\beta^{*}$ exists and is unique, margin condition (MC) is satisfied for $\varepsilon,\delta,\vartheta>0$, non-negative definite matrix $H\in R^{p_{n}\times p_{n}}$ and let

[TABLE]

where $\tilde{C}=2/(4+\varepsilon).$ Then:

[TABLE]

Proof.

Let:

[TABLE]

Lemmas 3 and 2 imply that:

[TABLE]

∎

Corollary 1**.**

(Separation property) If assumptions of Theorem 1 are satisfied, $\log p_{n}=o(\lambda^{2}n)$ and $\kappa_{H}(\varepsilon)>d$ for some $d,\varepsilon>0$ for large $n$, $|s^{*}|\lambda=o(\min\{\beta^{*}_{min},1\}),$ then

[TABLE]

Moreover

[TABLE]

Proof.

First part of the corollary follows directly from Theorem 1. Now we prove that condition $||\hat{\beta}_{L}-\beta^{*}||_{1}\leq\beta^{*}_{min}/2$ implies separation property

[TABLE]

Observe that for all $j\in\{1,\ldots,p_{n}\}$ we have:

[TABLE]

If $j\in s^{*},$ then using triangle inequality yields:

[TABLE]

Hence from the above inequality and (21) we obtain for $j\in s^{*}$:

[TABLE]

If $j\in s^{*c},$ then $\beta^{*}_{j}=0$ and (21) takes the form:

[TABLE]

This ends the proof. ∎

4 GIC consistency for a a general loss function and random prdictors

Theorems 2 and 3 state probability inequalities related to behaviour of GIC on supersets and on subsets of $s^{*}$, respectively. Corollaries 2 and 3 present asymptotic conditions for GIC consistency in the aforementioned situations. Corollary 4 gathers conclusions of Theorem 1 and Corollaries 1, 2 and 3 to show consistency of SS procedure (see Pokarowski and Mielniczuk (2015)) in case of subgaussian variables.

Theorem 2**.**

Assume that $\rho(\cdot,y)$ is convex, Lipschitz function with constant $L>0$, $X_{ij}\sim Subg(\sigma_{jn}^{2}),$ condition $C_{\epsilon}(w)$ holds for some $\epsilon,\theta>0$ and for every $w\subseteq\{1,\ldots,p_{n}\}$ such that $|w|\leq k_{n}$. Then for any $r<\epsilon$ we have:

[TABLE]

Proof.

If $s^{*}\subset w\in\mathcal{M}$ and $\hat{\beta}(w)-\beta^{*}\in B_{2}(r)$ then in view of inequalities $R_{n}(\hat{\beta}(s^{*}))\leq R_{n}(\beta^{*})$ and $R(\beta^{*})\leq R(b)$ we have:

[TABLE]

Note that $a_{n}(|w|-|s^{*}|)\geq a_{n}.$ Hence, if we have for some $w\supset s^{*}$: $GIC(w)\leq GIC(s^{*})$ then we obtain $nR_{n}(\hat{\beta}(s^{*}))-nR_{n}(\hat{\beta}(w)))\geq a_{n}(|w|-|s^{*}|)$ and from the above inequality we have $S_{1}(r)\geq{a_{n}}/{n}$. Furthermore, if $\hat{\beta}(w)-\beta^{*}\in B_{2}(r)^{c}$ and $r<\epsilon,$ then consider:

[TABLE]

where $u=r/(r+||\hat{\beta}(w)-\beta^{*}||_{2})$. Then

[TABLE]

as function $x/(x+r)$ is increasing with respect to $x$ for $x>0$. Moreover, we have $||v-\beta^{*}||_{2}\leq r<\epsilon$. Hence, in view of $C_{\epsilon}(w)$ condition we get:

[TABLE]

From convexity of $R_{n}$ we have:

[TABLE]

We observe that $\operatorname*{supp}v\subseteq\operatorname*{supp}\hat{\beta}(w)\cup\operatorname*{supp}\beta^{*}\subseteq w$, hence $v\in D_{1}$. Finally, we have:

[TABLE]

Hence we obtain the following sequence of inequalities:

[TABLE]

∎

Corollary 2**.**

Assume that $\rho(\cdot,y)$ is convex, Lipschitz function with constant $L>0$, $X_{ij}\sim Subg(\sigma_{jn}^{2}),$ condition $C_{\epsilon}(w)$ holds for some $\epsilon,\theta>0$ and for every $w\subseteq\{1,\ldots,p_{n}\}$ such that $|w|\leq k_{n}$, $k_{n}\ln(p_{n}\vee 2)=o(n)$ and $k_{n}\ln(p_{n}\vee 2)=o(a_{n})$. Then we have

[TABLE]

Proof.

We take: $r_{n}=C_{n}\sqrt{\frac{k_{n}\ln(p_{n}\vee 2)}{n}}$, where

[TABLE]

We observe that $C_{n}\rightarrow+\infty$, $r_{n}\leq\sqrt[4]{\frac{k_{n}\ln(p_{n}\vee 2)}{n}}\rightarrow 0$ and

[TABLE]

In view of Theorem 2 we have for sufficiently large $n$ such that $r_{n}<\epsilon$ holds:

[TABLE]

∎

The most restrictive condition of Corollary 2 is $k_{n}\ln(p_{n}\vee 2)=o(a_{n})$. We note that in the case when $p_{n}\geq n$ and $k_{n}=d$, EBIC penalty defined above corresponds to the borderline of this condition. Theorem 3 is an analogue of Theorem 2 for subsets of $s^{*}$.

Theorem 3**.**

Assume that $\rho(\cdot,y)$ is convex, Lipschitz function with constant $L>0$, $X_{ij}\sim Subg(\sigma_{jn}^{2}),$ condition $C_{\epsilon}(s^{*})$ holds for some $\epsilon,\theta>0$ and $8a_{n}|s^{*}|\leq\theta n\min\{\epsilon^{2},\beta^{*2}_{min}\}$. Then we have:

[TABLE]

Proof.

Suppose that for some $w\subset s^{*}$ we have $GIC(w)\leq GIC(s^{*})$. This is equivalent to:

[TABLE]

In view of inequalities $R_{n}(\hat{\beta}(s^{*}))\leq R_{n}(\beta^{*})$ and $a_{n}(|w|-|s^{*}|)\geq-a_{n}|s^{*}|$ we obtain:

[TABLE]

Let $v=u\hat{\beta}(w)+(1-u)\beta^{*}$ for some $u\in[0,1]$ to be specified later. From convexity of $\rho$ we consider:

[TABLE]

We consider two cases separately:

1. $\beta^{*}_{min}>\epsilon$.

First we observe

[TABLE]

what follows from our assumption. Let $u=\epsilon/(\epsilon+||\hat{\beta}(w)-\beta^{*}||_{2})$ and

[TABLE]

Note that $||\hat{\beta}(w)-\beta^{*}||_{2}\geq||\beta^{*}_{s^{*}\setminus w}||_{2}\geq\beta_{min}^{*}$. Then, as function $d(x)=x/(x+c)$ is increasing and bounded from above by $1$ for $x,c>0$, we obtain:

[TABLE]

Hence, in view of $C_{\epsilon}(s^{*})$ condition we have:

[TABLE]

Using (23)-(25) and above inequality yields:

[TABLE]

Thus, in view of Lemma 2, we obtain:

[TABLE]

2. $\beta^{*}_{min}\leq\epsilon$.

In this case we take $u=\beta^{*}_{min}/(\beta^{*}_{min}+||\hat{\beta}(w)-\beta^{*}||_{2})$ and define $v$ as in (25). Analogously as in (26), we have:

[TABLE]

Hence, in view of $C_{\epsilon}(s^{*})$ condition we have:

[TABLE]

Using (23) and above inequality yields:

[TABLE]

Thus, in view of Lemma 2, we obtain:

[TABLE]

By combining (27) and (28) the theorem follows. ∎

Corollary 3**.**

Assume that loss $\rho(\cdot,y)$ is convex, Lipschitz function with constant $L>0$, $X_{ij}\sim Subg(\sigma_{jn}^{2}),$ condition $C_{\epsilon}(s^{*})$ holds for some $\epsilon,\theta>0$ and $a_{n}|s^{*}|=o(n\min\{1,\beta^{*}_{min}\}^{2})$, then

[TABLE]

Proof.

First we observe as $a_{n}\to\infty$

[TABLE]

implies

[TABLE]

and thus in view of Theorem 3 we have

[TABLE]

∎

5 Selection consistency of SS procedure

In this section we combine the results of the two previous sections to establish consistency of a two-step SS procedure. It consists in construction of a nested family of models ${\cal M}$ using magnitude of Lasso coefficients and then finding the minimizer of GIC over this family. As ${\cal M}$ is data-dependent, in order to establish consistency of the procedure we use Corollaries 2 and 3 in which the minimizer of GIC is considered over all subsets and supersets of $s^{*}$.

SS (Screening and Selection) procedure is defined as follows:

Choose some $\lambda>0$. 2. 2.

Find $\hat{\beta}_{L}=\operatorname*{arg\min}\limits_{b\in R^{p_{n}}}R_{n}(b)+\lambda||b||_{1}$. 3. 3.

Find $\hat{s}_{L}=\operatorname*{supp}\hat{\beta}_{L}=\{j_{1},\ldots,j_{k}\}$ such that $|\hat{\beta}_{L,j_{1}}|\geq\ldots\geq|\hat{\beta}_{L,j_{k}}|>0$ and $j_{1},\ldots,j_{k}\in\{1,\ldots,p_{n}\}$. 4. 4.

Define $\mathcal{M}_{SS}=\{\emptyset,\{j_{1}\},\{j_{1},j_{2}\},\ldots,\{j_{1},j_{2},\ldots,j_{k}\}\}$. 5. 5.

Find $\hat{s}^{*}=\operatorname*{arg\min}\limits_{w\in\mathcal{M}_{SS}}GIC(w)$.

SS procedure is a modification of SOS procedure in Pokarowski et al. (2018) designed for GLMs for which ordering of variables is based on p-values of corresponding significance tests for them. Since additional ordering is omitted in the proposed modification we compressed the name to SS.

Corollary 4 and Remark 1 describe the situations when SS procedure is selection consistent. In it we use the assumptions imposed in Sections 2 and 3 together with an assumption that support of $s^{*}$ contains no more than $k_{n}$ elements, where $k_{n}$ is some deterministic sequence of integers.

Corollary 4**.**

Assume that $\rho(\cdot,y)$ is convex, Lipschitz function with constant $L>0$, $X_{ij}\sim Subg(\sigma_{jn}^{2})$ and $\beta^{*}$ exists and is unique. If $k_{n}\in N_{+}$ is some sequence, margin condition (MC) is satisfied for some $\vartheta,\delta,\varepsilon>0$, condition $C_{\epsilon}(w)$ holds for some $\epsilon,\theta>0$ and for every $w\subseteq\{1,\ldots,p_{n}\}$ such that $|w|\leq k_{n}$, $\mathcal{M}_{SS}$ is nested family constructed in the step 4 of SS procedure and the following conditions are fulfilled:

• •

$|s^{*}|\leq k_{n}$,

• •

$P(\forall w\in\mathcal{M}_{SS}:|w|\leq k_{n})\rightarrow 1$,

• •

$\liminf\limits_{n}\kappa_{H}(\varepsilon)>0$* for some $\varepsilon>0$, where $H$ is non-negative definite matrix and $\kappa_{H}(\varepsilon)$ is defined in (7),*

• •

$\log(p_{n})=o(n\lambda^{2})$,

• •

$k_{n}\lambda=o(\min\{\beta^{*}_{min},1\})$,

• •

$k_{n}\log p_{n}=o(n)$,

• •

$k_{n}\log p_{n}=o(a_{n})$,

• •

$a_{n}k_{n}=o(n\min\{\beta^{*}_{min},1\}^{2})$,

then for SS procedure we have

[TABLE]

Proof.

In view of Corollary 1 it follows from separation property (21) we obtain $P(s^{*}\in\mathcal{M}_{SS})\rightarrow 1$. Let:

[TABLE]

Then we have again from the fact that $A_{2}\cap B=\emptyset$, union inequality and Corollary 2:

[TABLE]

In the analogous way, using $|s^{*}|\leq k_{n}$ and Corollary 3 yields:

[TABLE]

Now, observe that in view of definition of $\hat{s}^{*}$ and union inequality:

[TABLE]

Thus $P(\hat{s}^{*}=s^{*})\rightarrow 1$ in view of above inequality, (29) and (30). ∎

Consider now the case of semiparametric model defined in (4). Then it is known (cf Brillinger (1982), Ruud (1983) and Li and Duan (1989)) that provided $X$ has regressions satisfying

[TABLE]

where $\beta$ is the true parameter, then $\beta^{*}=\eta\beta$ and $\eta\neq 0$ if ${\rm Cov}(Y,X)\neq 0$. The linear regressions condition (31) is satisfied e.g. by eliptically contoured distribution, in particular by multivariate normal. We refer also to Kubkowski and Mielniczuk (2017)) for discussion and up-to date references to this problem. In the discussed case $s^{*}=s$ and we can state the following result.

Corollary 5**.**

Assume that (31) and assumption of Corollary 4 are satisfied. Then $P(\hat{s}^{*}=s)\to 1$.

Remark 1**.**

If $p_{n}=O(e^{cn^{\gamma}})$ for some $c>0$, $\gamma\in(0,1/2)$, $\xi\in(0,0.5-\gamma)$, $u\in(0,0.5-\gamma-\xi)$, $k_{n}=O(n^{\xi})$, $\lambda=C_{n}\sqrt{\log(p_{n})/n}$, $C_{n}=O(n^{u})$, $C_{n}\rightarrow+\infty$, $n^{-\frac{\gamma}{2}}=O(\beta^{*}_{min})$, $a_{n}=dn^{\frac{1}{2}-u}$, then assumptions imposed on asymptotic behaviour of parameters in Corollary 4 are satisfied.

Remark 2**.**

We note that in order to apply Corollary 4 to two-step procedure based on Lasso it is required that $|s^{*}|\leq k_{n}$ and that the support of Lasso estimator with probability tending to $1$ contains no more than $k_{n}$ elements. Some results bounding $|\operatorname*{supp}\hat{\beta}_{L}|$ are available for deterministic $X$ (see Huang et al. (2008)) and for random $X$ (see Tibshirani (2013)), but they are too weak to be useful for EBIC penalties. The other possibility to prove consistency of two-step procedure is to modify it in the first step by using thresholded Lasso (see Zhou (2010)) corresponding to $k_{n}^{\prime}$ largest Lasso coefficients where $k_{n}^{\prime}\in N$ is such that $k_{n}=o(k_{n}^{\prime})$.

6 Numerical experiments

6.1 Selection procedures

In performed simulations we have implemented modifications of SS procedure introduced in Section 5, as the original procedure is defined for a single $\lambda$ only. In practice it is generally more convenient to consider some sequence $\lambda_{1}>\ldots>\lambda_{m}>0$ instead of only one $\lambda$ in the first step in order to avoid choosing ’the best’ $\lambda$. For the chosen sequence $\lambda_{1},\ldots,\lambda_{m}$, we construct corresponding families $\mathcal{M}_{1},\ldots,\mathcal{M}_{m}$ analogously to $\mathcal{M}$ in the step 4 of SS procedure. Thus we arrive here at the following SSnet procedure, which is the modification of SOSnet procedure in Pokarowski et al. (2018). Below $\tilde{b}$ is a vector $b$ with first coordinate corresponding to intercept omitted, $b=(b_{0},{\tilde{b}}^{T})^{T})$:

Choose some $\lambda_{1}>\ldots>\lambda_{m}>0$. 2. 2.

Find $\hat{\beta}_{L}^{(i)}=\operatorname*{arg\min}\limits_{b\in R^{p_{n}+1}}R_{n}(b)+\lambda_{i}||\tilde{b}||_{1}$ for $i=1,\ldots,m$. 3. 3.

Find $\hat{s}_{L}^{(i)}=\operatorname*{supp}\hat{\tilde{\beta}}_{L}^{(i)}=\{j_{1}^{(i)},\ldots,j_{k_{i}}^{(i)}\}$ where $j_{1}^{(i)},\ldots,j_{k_{i}}^{(i)}$ are such that $|\hat{\beta}_{L,j_{1}^{(i)}}^{(i)}|\geq\ldots\geq|\hat{\beta}_{L,j_{k}^{(i)}}^{(i)}|>0$ for $i=1,\ldots,m$. 4. 4.

Define $\mathcal{M}_{i}=\{\{j_{1}^{(i)}\},\{j_{1}^{(i)},j_{2}^{(i)}\},\ldots,\{j_{1}^{(i)},j_{2}^{(i)},\ldots,j_{k_{i}}^{(i)}\}\}$ for $i=1,\ldots,m$. 5. 5.

Define $\mathcal{M}=\{\emptyset\}\cup\bigcup\limits_{i=1}^{m}\mathcal{M}_{i}$. 6. 6.

Find $\hat{s}^{*}=\operatorname*{arg\min}\limits_{w\in\mathcal{M}}GIC(w)$, where

[TABLE]

Instead of constructing families $\mathcal{M}_{i}$ for each $\lambda_{i}$ in SSnet procedure, $\lambda$ can be chosen by cross-validation using 1SE rule (see Friedman et al. (2010)) and then proceed as in SS procedure. We call this procedure SSCV.

The last procedure considered has been introduced by Fan and Tang in Fan and Tang (2013) and is Lasso procedure with penalty parameter $\hat{\lambda}$ chosen in a data-dependent way as for SSCV. Namely, it is the minimizer of GIC criterion with $a_{n}=\log(\log n)\cdot\log p_{n}$ for which ML estimator has been replaced by Lasso estimator with penalty $\lambda$. Once $\hat{\beta}_{L}(\hat{\lambda}_{L})$ is calculated then $\hat{s}^{*}$ is defined as its support. The procedure is called LFT in the sequel.

We list below versions of the above procedures along with R packages, which were used to choose sequence $\lambda_{1},\ldots,\lambda_{m}$ and computation of Lasso estimator. The following packages were chosen based on selection performance after initial tests for each loss and procedure:

• •

SSnet with logistic or quadratic loss: ncvreg,

• •

SSCV or LFT with logistic or quadratic loss: glmnet,

• •

SSnet, SSCV or LFT with Huber loss (cf Yi and Huang (2017)): hqreg.

The following functions which were used to optimize $R_{n}$ in GIC minimization step for each loss:

• •

logistic loss: glm.fit (package stats),

• •

quadratic loss: .lm.fit (package stats),

• •

Huber loss: rlm (package rlm).

Before applying investigated procedures, each column of matrix

[TABLE]

was standardized as $\hat{\beta}_{L}$ depends on scaling of predictors. We set length of $\lambda_{i}$ sequence to $m=20$. Moreover, in all procedures we considered only $\lambda_{i}$ for which $|\hat{s}_{L}^{(i)}|\leq n$. It is due to the fact that when $|\hat{s}_{L}^{(i)}|>n$ Lasso and ML solutions are not unique (see Rosset et al. (2004), Tibshirani (2013)). For Huber loss we set parameter $\delta=1/10$ (see Yi and Huang (2017)). Number of folds in SSCV was set to $K=10$.

Each simulation run consisted of $L$ repetitions, during which samples $\mathbb{X}_{k}=(X_{1}^{(k)},\ldots,X_{n}^{(k)})^{T}$ and $\mathbf{Y}_{k}=(Y_{1}^{(k)},\ldots,Y_{n}^{(k)})^{T}$ were generated for $k=1,\ldots,L$. For $k$-th sample $(\mathbb{X}_{k},\mathbf{Y}_{k})$ estimator $\hat{s}_{k}^{*}$ of set of active predictors is obtained by a given procedure as the support of $\hat{\tilde{\beta}}(\hat{s}_{k}^{*})$, where

[TABLE]

is ML estimator for $k$-th sample. We denote by $\mathcal{M}^{(k)}$ is the family $\mathcal{M}$ obtained by a given procedure for $k$-th sample.

In our numerical experiments we have computed the following measures of selection performance which gauge co-direction of true parameter $\beta$ and $\hat{\beta}$ and the interplay between $s^{*}$ and $\hat{s}^{*}$:

• •

$ANGLE=\frac{1}{L}\sum\limits_{k=1}^{L}\arccos|\cos\angle(\tilde{\beta}_{0},\hat{\tilde{\beta}}(\hat{s}_{k}^{*}))|$, where

[TABLE]

and we let $\cos\angle(\tilde{\beta},\hat{\tilde{\beta}}(\hat{s}_{k}^{*}))=0$, if $||\tilde{\beta}||_{2}||\hat{\tilde{\beta}}(\hat{s}_{k}^{*})||_{2}=0$,

• •

$P_{inc}=\frac{1}{L}\sum\limits_{k=1}^{L}I(s^{*}\in\mathcal{M}^{(k)})$,

• •

$P_{equal}=\frac{1}{L}\sum\limits_{k=1}^{L}I(\hat{s}_{k}^{*}=s^{*})$.

• •

$P_{supset}=\frac{1}{L}\sum\limits_{k=1}^{L}I(\hat{s}_{k}^{*}\supseteq s^{*})$.

6.2 Regression models considered

6.2.1 Model M1

In order to investigate behaviour of two-step procedure under misspecification we considered two similar models with different sets of predictors. As sets of predictors differ, this results in correct specification of the first model (model M1) and misspecification of the second (Model M2).

Namely, we generated $n$ observations $(X_{i},Y_{i})\in R^{p+1}\times\{0,1\}$ for $i=1,\ldots,n$ such that:

[TABLE]

where $Z_{i}=(Z_{i1},\ldots,Z_{ip})^{T}\sim\mathcal{N}_{p}(0_{p},\Sigma)$, $\Sigma=[\rho^{|i-j|}]_{i,j=1,\ldots,p}$ and $\rho\in(-1,1)$. We consider response function $q(x)=q_{L}(x^{3})$ for $x\in R$, $s=\{1,2\}$ and $\beta_{s}=(1,1)^{T}$. This means that:

[TABLE]

We observe that the above binary model is well specified with respect to family of fitted logistic models. Hence $s^{*}=\{6,7,8,9\}$ and $\beta^{*}_{s^{*}}=(3,3,1,1)^{T}$ are respectively set of active predictors and non-zero coefficients of projection onto family of logistic models.

We considered the following parameters in numerical experiments: $n=500,p=150,\rho\in\{-0.9+0.15\cdot k\colon\leavevmode\nobreak\ k=0,1,\ldots,12\}$ and $L=500$ - number of generated data sets for each combination of parameters. We investigated procedures SSnet, SSCV and LFT using logistic, quadratic and Huber (cf Yi and Huang (2017)) loss functions. For procedures SSnet and SSCV we used GIC penalties with:

• •

$a_{n}=\log n$ (BIC),

• •

$a_{n}=\log n+2\log p_{n}$ (EBIC1).

6.2.2 Model M2

We generated $n$ observations $(X_{i},Y_{i})\in R^{p+1}\times\{0,1\}$ for $i=1,\ldots,n$ such that $X_{i}=(X_{i0},X_{i1},\ldots,X_{ip})^{T}$ and $(X_{i1},\ldots,X_{ip})^{T}\sim\mathcal{N}_{p}(0_{p},\Sigma)$, $\Sigma=[\rho^{|i-j|}]_{i,j=1,\ldots,p}$ and $\rho\in(-1,1)$. Response function is $q(x)=q_{L}(x^{3})$ for $x\in R$, $s=\{1,2\}$ and $\beta_{s}=(1,1)^{T}$. This means that:

[TABLE]

This model in comparison to the one presented in Section 6.2.1 does not contain monomials of $X_{i1}$ and $X_{i2}$ of degree higher than $1$ in its set of predictors. We observe that this binary model is missspecified with respect to fitted family of logistic models, because $q(x_{i1}+x_{i2})\not\equiv q_{L}(\beta^{T}x_{i})$ for any $\beta\in R^{p+1}$. However, in this case linear regressions condition (31) is satisfied for $X$, as it follows normal distribution (see Kubkowski and Mielniczuk (2018),Li and Duan (1989)) . Hence in view of Proposition 3.8 in Kubkowski and Mielniczuk (2017) we have $s^{*}_{log}=\{1,2\}$ and $\beta^{*}_{log,s^{*}_{log}}=\eta(1,1)^{T}$ for some $\eta>0$. Parameters $n,p,\rho$ as well as $L$ were chosen as for model M1.

6.2.3 Results for models M1 and M2

We look first at behaviour of $P_{inc}$, $P_{equal}$ and $P_{supset}$ for the considered procedures. We observe that values of $P_{inc}$ for SSCV and SSnet are close to $1$ for low correlations in model M2 for every tested loss (see Figure 1). In model M1 $P_{inc}$ attains the largest values for SSnet procedure and logistic loss for low correlations - this is due to the fact that in most cases corresponding family $\mathcal{M}$ is the largest among the families created by considered procedures. $P_{inc}$ is close to [math] in model M1 for quadratic and Huber loss, what results in low values of the remaining indices. This may be due to strong dependences between predictors in model M1, note e.g. that we have $\operatorname*{Cor}(X_{i1},X_{i8})=3/\sqrt{15}\approx 0.77$. It is seen that in model M1 inclusion probability $P_{inc}$ is much lower than in model M2 (except for negative correlations). It it also seen that $P_{inc}$ for SSCV is larger than for LFT and LFT fails with respect to $P_{inc}$ in M1.

In model M1 the largest values $P_{equal}$ are attained for SSnet with BIC penalty, then for SSCV with EBIC1 penalty (see Figure 2). In the model M2 $P_{equal}$ is close to $1$ for SSnet and SSCV with EBIC1 penalty and was much larger than $P_{equal}$ for the corresponding versions using BIC penalty. We also note that choice of loss was relevant only for larger correlations. These results confirm theoretical result of Theorem 2.1 in Li and Duan (1989) which show that collinearity holds for broad class of loss function. We observe also that although in the model M2 remaining procedures do not select $s^{*}$ with high probability, they select its superset, what is indicated by values of $P_{supset}$ (see Figure 3). This analysis is confirmed by an analysis of $ANGLE$ measure (see Figure 4), which attains values close to [math], when $P_{supset}$ is close to $1$. Low values of $ANGLE$ measure mean that estimated vector $\hat{\tilde{\beta}}(\hat{s}_{k}^{*})$ is approximately proportional to $\tilde{\beta}$, what was the case for M2 model, where normal predictors satisfy linear regressions condition. Note that the angles of $\hat{\tilde{\beta}}(\hat{s}_{k}^{*})$ and $\tilde{\beta}^{*}$ in M1 significantly differ despite the fact that M1 is well specified. Also, for the best performing procedures in both models and any loss considered, $P_{equal}$ was much larger in M2 than in M1, despite the fact that the latter is correctly specified. This shows that choosing a simple misspecified model which retains crucial characteristics of the well specified large model instead of the latter might be beneficial.

In model M1 procedures with BIC penalty performed better than those with EBIC1 penalty, however the gain for $P_{equal}$ was much smaller than the gain when using EBIC1 in M2. LFT procedure performed poorly in model M1 and reasonably well in model M2. The overall winner in both models is SSnet. SSCV performs only slightly worse than SSnet in M2 but performs significantly worse in M1.

Analysis of computing times of 1st and 2nd stage of each procedure shows that SSnet procedure creates large families $\mathcal{M}$ and GIC minimization becomes computationally intensive. We also observe that the first stage for SSCV is more time consuming than for SSnet, what is caused by multiple fitting of Lasso in cross-validation. However, SSCV is much faster than SSnet in the second stage.

We conclude that in the considered experiments SSnet with EBIC1 penalty works the best in most cases, however even for the winning procedure strong dependence of predictors results in deterioration of its performance. It is also clear from our experiments that a choice of GIC penalty is crucial for its performance. Modification of SS procedure which would perform satisfactorily for large correlations is still an open problem.

7 Appendix

Proof of Lemma 1:

Proof.

Firstly, observe that function $R_{n}$ is convex as $\rho$ is convex. Moreover, from the definition of $\hat{\beta}_{L}$ we get the inequality:

[TABLE]

We note that $v-\beta^{*}\in B_{1}(r),$ as we have:

[TABLE]

By definition of $W_{n},$ convexity of $R_{n}$, (33) and definition of $S$ we have:

[TABLE]

From the convexity of $l_{1}$ norm, (34), (32), $||\beta^{*}||_{1}=||\beta^{*}_{s^{*}}||_{1}$ and triangle inequality it follows that:

[TABLE]

Hence:

[TABLE]

∎

Lemma 4**.**

Assume that $S\sim Subg(\sigma^{2})$ and $T$ is random variable such that $|T|\leq M,$ where $M$ is some positive constant and $S$ and $T$ are independent. Then $ST\sim Subg(M^{2}\sigma^{2}).$

Proof.

Observe that:

[TABLE]

∎