Selecting Penalty Parameters of High-Dimensional M-Estimators using Bootstrapping after Cross-Validation

Selecting Penalty Parameters of High-Dimensional M-Estimators using Bootstrapping after Cross-Validation

Keywords: Penalty parameter selection, penalized M-estimation, high-dimensional models, sparsity, cross-validation, bootstrap, inference, one-step debiasing.

Introduction

High-dimensional models have attracted substantial attention both in the econometrics and in the statistics/machine learning literature, see e.g. belloni2018highdimensional and hastie_statistical_2015, and $\ell_{1}$-penalized estimators have emerged among the most useful methods for learning parameters of such models. However, implementing these estimators requires a choice of the penalty parameter and with few notable exceptions, e.g. $\ell_{1}$-penalized linear mean, quantile and logit regression estimators, the choice of this penalty parameter in practice often remains unclear. In this paper, we develop a new method to choose the penalty parameter in the context of $\ell_1$-penalized M-estimation and show that our method leads to precise estimation and inference in a large variety of models.

We consider a model where the true value $\boldsymbol{\theta}_{0}$ of some parameter $\boldsymbol{\theta}$ is given by the solution to an optimization problem

where $m:\mathbb{R}\times \mathcal Y \to \mathbb{R}$ is a known (potentially non-smooth) loss function that is convex in its first argument, $\boldsymbol{X} = (X_{1},\dots,X_{p})^{\top}\in\mathcal{X}\subseteq\mathbb{R}^{p}$ a vector of candidate regressors, $\boldsymbol{Y}\in\mathcal{Y}$ one or more outcome variables, and $\Theta\subseteq\mathbb{R}^{p}$ a convex parameter space. Prototypical loss functions are square loss and negative log-likelihood, but the framework ((ref)) also covers many other cross-sectional models and associated modern as well as classical estimation approaches including logit and probit models, logistic calibration tan2020regularized, covariate balancing imai_covariate_2014, and expectile regression newey_asymmetric_1987. It also subsumes approaches to estimation of panel-data models such as the fixed-effects/conditional logit for binary outcomes rasch1960probabilistic, trimmed least-absolute-deviations and trimmed least-squares for censored outcomes honore_trimmed_1992, and partial likelihood approaches to heterogeneous panel models for duration chamberlain_heterogeneity_1985. We detail some of these examples in Section (ref).

For the purpose of estimation, we assume access to a sample $\{(\boldsymbol{X}_{i},\boldsymbol{Y}_{i})\}_{i=1}^n$ of $n$ independent observations from the distribution $P$ of the pair $(\boldsymbol{X},\boldsymbol{Y})$, where the number $p$ of candidate regressors in each $\boldsymbol{X}_{i}=(X_{i,1},\dots,X_{i,p})^{\top}$ may be (potentially much) larger than the sample size $n$, meaning that we cover high-dimensional models. Following the literature on high-dimensional models, we assume that the vector $\boldsymbol{\theta}_{0}=(\theta_{0,1},\dots,\theta_{0,p})^{\top}$ is at least approximately (also known as “weakly”) sparse. While we postpone a formal definition to Section (ref), approximate sparsity captures the idea that, even though the number of candidate regressors $p$ can be very large, the number of relevant regressors may be substantially smaller. In the simplest case, known as exact (or “strong”) sparsity, this assumption amounts to the number of non-zeros in $\boldsymbol{\theta}_0$ being much smaller than $n$. Approximate sparsity relaxes this idea to allow possibly many---but typically small---non-zeros. With sparsity in mind, we study the sparsity encouraging $\ell_{1}$-penalized M-estimator ($\ell_{1}$-ME)

where $\Vert\boldsymbol{\theta}\Vert_{1}=\sum_{j=1}^{p}\vert\theta_{j}\vert$ denotes the $\ell_{1}$ norm of $\boldsymbol{\theta}$, and $\lambda\in[0,\infty)$ is a penalty parameter.\footnote{Throughout the main text, we implicitly assume that an estimator exists. Simple conditions under which $\widehat{\Theta}(\lambda)$ is non-empty (and related properties) are given in Appendix (ref).} We also study the post-$\ell_1$-penalized M-estimator (post-$\ell_1$-ME), which refits the coefficients of the variables selected by $\ell_1$-ME without the penalty in the criterion function in (ref).

Implementing the estimator $\widehat{\boldsymbol{\theta}}(\lambda)$ requires choosing $\lambda$. To do so, we first extend a probabilistic bound from belloni_l1-penalized_2011, obtained for $\ell_{1}$-penalized quantile regression, to our general $\ell_{1}$-penalized M-estimation setting ((ref)). (See also negahban_unified_2012 for independently developed and closely related results.) The bound, which we state in Section (ref), yields a general principle to choose $\lambda$. In particular, it suggests that, for an arbitrary choice of $c_{0}\in(1,\infty)$, one should choose $\lambda$ as small as possible subject to the constraint that the event

occurs with probability approaching one, where $m_1'$ denotes the partial derivative of the loss function with respect to its first argument. We therefore wish to set $\lambda=c_{0}q_n(1-\alpha)$, where

for some small user-specified probability tolerance level $\alpha = \alpha_n \to 0$ as $n\to\infty$. This choice, however, is typically infeasible since the random variable in ((ref)) depends on the unknown $\boldsymbol{\theta}_{0}$. We thus have a vicious circle: to choose $\lambda$, we need an estimator of $\boldsymbol{\theta}_{0}$, but to estimate $\boldsymbol{\theta}_{0}$, we need to choose $\lambda$. In this paper, we offer a solution to this problem, which constitutes our key contribution.

To obtain our solution, we show that even though (as we discuss below) the estimator $\widehat{\boldsymbol{\theta}}(\lambda)$ based on $\lambda$ chosen by cross-validation or its variants is generally difficult to analyze, it can be used to construct provably good, in a certain sense, estimators of the random vectors $m_{1}'(\boldsymbol{X}_{i}^{\top}\boldsymbol{\theta}_{0},\boldsymbol{Y}_{i})\boldsymbol{X}_{i}$. We are then able to derive an estimator, say $\widehat{q}(1-\alpha)$, of $q_n(1-\alpha)$ via bootstrapping, as discussed in belloni2018highdimensional, and to set $\lambda=c_{0}\widehat{q}(1-\alpha)$, which we refer to as the bootstrap-after-cross-validation (BCV) method to choose $\lambda$. This method is computationally rather straightforward, applicable in a wide variety of models, and non-conservative in the sense that it gives $\lambda$ such that $\lambda\approx c_{0}q_n(1-\alpha)$ rather than $\lambda \gg c_0q_n(1-\alpha)$. We derive convergence rates of $\ell_1$-ME and post-$\ell_1$-ME based on this choice of $\lambda$ in Section (ref). In addition, we show in Section (ref) that, upon debiasing via the double machine learning approach, these estimators yield simple inference procedures.

The main alternatives to our method are cross-validation and related sample-splitting techniques. One of the main complications with these methods is that they are difficult to analyze, at least in some important dimensions. Sample-splitting techniques yield bounds on the $\ell_{2}$ estimation error $\|\widehat{\boldsymbol{\theta}}(\lambda)-\boldsymbol{\theta}_{0}\|_{2}$, see e.g. lecue_mitchell, but not on the $\ell_{1}$ estimation error $\|\widehat{\boldsymbol{\theta}}(\lambda)-\boldsymbol{\theta}_{0}\|_{1}$.\footnote{Any two norms on a fixed and finite-dimensional space are equivalent. However, the equivalence constants generally depend on the dimension (here $p$), which makes translation of error bounds for one norm into another a non-trivial manner when the dimension is growing.} In contrast, our method gives bounds on both $\ell_2$ and $\ell_1$ estimation errors. An $\ell_1$ error bound is crucial when we are interested in estimating dense functionals $\boldsymbol{a}^{\top}\boldsymbol{\theta}_{0}$ of $\boldsymbol{\theta}_0$ with $\boldsymbol{a}\in\mathbb{R}^{p}$ being a vector of loadings with many non-zero components; see belloni2018highdimensional for details.\footnote{Dense functionals $\boldsymbol{a}^{\top}\boldsymbol{\theta}_0$ may appear in the analysis, for example, when the vector $\boldsymbol{X}$ consists of many dummy variables and we are interested in making comparisons between two cells, $(\boldsymbol{x}_2 - \boldsymbol{x}_1)^{\top}\boldsymbol{\theta}_0$, where $\boldsymbol{x}_1$ and $\boldsymbol{x}_2$ represent the first and the second cell, respectively. In such examples, we can guarantee that $(\boldsymbol{x}_2 - \boldsymbol{x}_1)^{\top}\widehat\boldsymbol{\theta}(\lambda)$ is close to $(\boldsymbol{x}_2 - \boldsymbol{x}_1)^{\top}\boldsymbol{\theta}_0$ only when $\|\widehat\boldsymbol{\theta}(\lambda) - \boldsymbol{\theta}_0\|_1$ is small.} Moreover, $\ell_1$ estimation error bounds are needed to perform inference on components of $\boldsymbol{\theta}_{0}$ as in Section (ref).\footnote{It is possible to replace the requirement on $\ell_1$ estimation error by the requirement on $\ell_2$ estimation error via cross-fitting, as in chernozhukov2018double. However, the combination of sample-splitting and cross-fitting would require splitting the original sample into at least three subsamples, which may not lead to accurate inference in moderate samples.} When $\lambda$ is selected by cross-validation, $\ell_{1}$ and $\ell_{2}$ estimation error bounds are typically both unknown. The only exception we are aware of is the linear mean regression model estimated by the LASSO. For this special case, bounds have been derived in chetverikov_cross-validated_2016 and miolane_montanari, but the bounds appearing in those references are less sharp than those provided here. Moreover, and crucially, cross-validation may lead to rather poor inference results, in the sense of bad size control, even in relatively large samples, and does not dominate our method even in terms of estimation errors; see our simulation results in Section (ref) for details.

Another alternative to our method is to base the penalty parameter choice on self-normalized moderate deviation (SNMD) theory, as proposed in belloni_sparse_2012 for the linear mean regression model and extended in belloni2016post to the logit model. This method is slightly conservative, in the sense that it gives $\lambda$ somewhat larger than $c_0 q_n(1-\alpha)$, but yields estimation and inference results that are comparable in quality with those produced by the BCV method. The SNMD method can be further extended to cover any Lipschitz-continuous loss function, but it is not clear how to extend it to a non-Lipschitz setting. For example, the SNMD method can be applied to the logit model but not to the probit model. In contrast, our BCV method is nearly universally applicable, and does not require Lipschitz continuity. We provide several other important examples where the loss function is not Lipschitz-continuous in Section (ref).

To showcase our method using real data, in Section (ref) we revisit the setting of fryer_jr_empirical_2019, who investigated racial differences in police use of force. We extend Fryer's regression analysis in two ways. First, we change the model from a binary logit to a binary probit, keeping the regressors as in Fryer's analysis, a relatively small list. Second, we add a large number of additional (technical) regressors resulting from interactions between the original regressors. The first change leads to a non-Lipschitz loss (the negative probit likelihood). The second change brings us into high-dimensional territory, causing classical methods to break down. Unlike existing methods, the methods developed in this paper can accommodate both challenges. Our analysis supports the conclusions of fryer_jr_empirical_2019 in showing that they are robust to model specification and a much larger set of candidate controls than originally considered.

The literature on learning parameters of high-dimensional models via $\ell_{1}$-penalized M-estimation is large. Instead of listing all existing papers, we therefore refer the interested reader to the excellent textbook treatment in wainwright_high-dimensional_2019 and focus here on only a few key references. van_de_geer_high-dimensional_2008,van_de_geer_estimation_2016 derives bounds on the estimation errors of general $\ell_{1}$-penalized M-estimators ((ref)) and provides some choices of the penalty parameter $\lambda$. However, her penalty formulas give values of $\lambda$ that are so large that the resulting estimators are typically trivial in moderate samples, with all coefficients being exactly zero. Recognizing this issue, van_de_geer_high-dimensional_2008 remarks that her results should only be seen as an indication that her theory has something to say about finite sample sizes, and that other methods to choose $\lambda$ should be used in practice. negahban_unified_2012 develop error guarantees in a very general setting, and when specialized to our setting ((ref)), their results become quite similar to those in our Theorem (ref). The same authors also note that a challenge to using these results in practice is that the random variable in ((ref)) is usually impossible to compute because it depends on the unknown vector $\boldsymbol{\theta}_{0}$ (ibid., p. 547). It is exactly this challenge that we overcome in this paper. belloni_l1-penalized_2011 study the high-dimensional quantile regression model and note that the distribution of the random variable in ((ref)) is in this case pivotal, making the choice of the penalty parameter simple. Similarly, wang_tuning-free_2020 study the high-dimensional mean regression model and show that one can obtain pivotality by replacing the square-loss function by Jaeckel's dispersion function, again making the choice of the penalty parameter simple. However, these are the only two settings we are aware of in which the distribution of the random variable in ((ref)) is pivotal.\footnote{With a known censoring propensity, the buchinsky1998alternative linear programming estimator for censored quantile regression boils down to a variant of quantile regression, thus leading to pivotality.} Finally, ninomiya_kawano consider information criteria for the choice of the penalty parameter $\lambda$ but focus on fixed-$p$ asymptotics, thus precluding high-dimensional models.

The rest of the paper is organized as follows. In Section (ref) we provide a portfolio of examples that constitute possible applications of our method. In Section (ref) we develop bounds on the estimation error of the $\ell_1$-ME, which motivate our method for choosing the penalty parameter. In Section (ref), we introduce the BCV penalty method and derive convergence rates for the resulting $\ell_1$-ME and post-$\ell_1$-ME. In Section (ref), we show how to perform inference on individual components of $\boldsymbol{\theta}_0$ via debiasing. In Section (ref), we present a simulation study shedding light on the finite-sample properties of our method and contrast it with cross-validation. Finally, in Section (ref), we apply our method to the empirical setting of fryer_jr_empirical_2019. All proofs are relegated to the Online Appendices.

Notation

The distribution $P$ of the pair $(\boldsymbol{X},\boldsymbol{Y})$ and features thereof, including the dimension $p$ of the vector $\boldsymbol{X}$, may change with the sample size $n$ (that is, we consider triangular array sampling and asymptotics), but we suppress this potential dependence whenever this does not cause confusion in order to simplify notation. We use $\mathrm{E}[f(\boldsymbol{X},\boldsymbol{Y})]$ (or $\mathrm{E}_{\boldsymbol{X},\boldsymbol{Y}}[f(\boldsymbol{X},\boldsymbol{Y})]$) to denote the expectation of a function $f$ of the pair $(\boldsymbol{X},\boldsymbol{Y})$ computed with respect to $P$, and we use $\mathbb{E}_{n}[f(\boldsymbol{X}_i,\boldsymbol{Y}_i)]:=n^{-1}\sum_{i=1}^{n}f(\boldsymbol{X}_i,\boldsymbol{Y}_i)$ to abbreviate the sample average. We use $\mathbb{R}$ and $\mathbb{N}$ to denote all real numbers and all positive integers $\{1,2,\dotsc\},$ respectively. For $k\in\mathbb{N}$, we write $\left[k\right]:=\{1,\dotsc,k\}$ for all positive integers up to and including $k$. When only a non-empty subset $I\subsetneq[n]$ is in use, we write $\mathbb{E}_{I}[f(\boldsymbol{X}_{i},\boldsymbol{Y}_i)]:=|I|^{-1}\sum_{i\in I}f(\boldsymbol{X}_{i},\boldsymbol{Y}_i)$ for the subsample average. For a set of indices $I\subseteq\left[n\right]$, we use $I^{c}$ to denote the elements of $\left[n\right]$ not in $I.$ For $k\in\mathbb{N}$, we use $\mathbf 0_k$ to denote the vector in $\mathbb{R}^k$ whose components are all zero. Given a vector $\boldsymbol{\delta}\in\mathbb{R}^{k},$ we denote its $\ell_{r}$ norms, $r\in [1,\infty]$, by $\left\lVert\boldsymbol{\delta}\right\rVert_{r}$. We write $\mathrm{supp}(\boldsymbol{\delta}):=\{j\in[k];\boldsymbol{\delta}_j\neq0\}$ for the support of $\boldsymbol{\delta}$, and use the $\ell_0$ “norm” $\|\boldsymbol{\delta}\|_0:=|\mathrm{supp}(\boldsymbol{\delta})|$ to denote the number of non-zero elements of $\boldsymbol{\delta}$, where $|J|$ denotes the cardinality of the set $J$. For any function $f\colon \mathbb{R}\times \mathcal Z\to\mathbb{R}$, whose first argument is a scalar, we use $f'_1$, $f''_{11}$ and $f'''_{111}$ to denote its partial derivatives with respect to the first argument of the first, second and third order, respectively. We abbreviate $a\lor b:=\max\{a,b\} $ and $a\wedge b:=\min\{ a,b\} $. Unless explicitly stated otherwise, limits are understood as $n\to\infty$. For numbers $a_n$ and positive numbers $b_n,n\in\mathbb{N},$ we write $a_n=o(1)$ if $a_n\to0$, and $a_n\lesssim b_n,$ if the sequence $a_n/b_n$ is bounded. For random variables $V_n$ and positive numbers $b_n,$ we write $V_n\lesssim_{\mathrm{P}} b_n,$ if the sequence $V_n/b_n$ is bounded in probability. We denote $\eta_n:=\sqrt{\ln(pn)/n}$. We use the word “constant” to refer to non-random quantities that do not depend on $n$. Finally, we take $n\geqslant3$ and $p\geqslant2$ throughout and introduce more notation as needed in the appendices.

Examples

In this section, we discuss a variety of models that fit into the M-estimation framework ((ref)) with the loss function $m(t,\boldsymbol{y})$ being convex in its first argument. The following examples cover both discrete and continuous outcomes in likelihood and non-likelihood settings with smooth as well as kinked loss functions. Additional examples can be found in Appendix (ref).

Non-Asymptotic Bounds on Estimation Error

In this section, we derive probabilistic bounds on the error of the $\ell_{1}$-ME ((ref)) in the $\ell_{1}$ and $\ell_{2}$ norms. The bounds reveal which quantities one needs to control in order to ensure good behavior of the estimator, motivating the choice of the penalty parameter $\lambda$ in the next section.

Our bounds will be based on the following assumptions. Since Assumptions (ref), (ref), and (ref) stated below are high level, we verify these assumptions under more low-level conditions in the familiar case of the linear model with square loss in Appendix (ref) and for all examples in Section (ref) in Appendix (ref).

Assumption (ref) is a minor regularity condition. Both convexity and interiority follow trivially in the case of a full parameter space $\Theta=\mathbb{R}^p$. Assumption (ref) is satisfied in all examples from the previous section, as discussed there. In the same examples, Assumption (ref) imposes minor integrability conditions on the random vectors $\boldsymbol{X}$ and $\boldsymbol{Y}$. In addition, in the case of Example (ref) with trimmed LAD loss function, this assumption requires that the conditional distribution of $Y_1 - Y_2$ given $(\boldsymbol{X}, Y_1 + Y_2 > 0)$ is continuous; see Appendix (ref) for details.

Further, define the excess risk function $\mathcal{E}:\Theta\to[0,\infty)$ by \[ \mathcal{E}\left(\boldsymbol{\theta}\right):=\mathrm{E}\left[m\left(\boldsymbol{X}^{\top}\boldsymbol{\theta},\boldsymbol{Y}\right)-m\left(\boldsymbol{X}^{\top}\boldsymbol{\theta}_{0},\boldsymbol{Y}\right)\right],\quad\boldsymbol{\theta}\in\Theta. \] By definition of $\boldsymbol{\theta}_0$ in (ref), this function is non-negative and takes value zero at $\boldsymbol{\theta} = \boldsymbol{\theta}_0$. The next assumption requires that it grows sufficiently fast as $\boldsymbol{\theta}$ moves away from $\boldsymbol{\theta}_0$.

In addition to some technical regularity conditions, this assumption requires the matrix $\mathrm{E}[\boldsymbol{X}\boldsymbol{X}^\top]$ to be non-singular, which means that there should be no perfect regressor multicollinearity in the population. In the context of Example (ref), it also requires $Y_1$ and $Y_2$ to be different with positive probability. Also, our formal analysis reveals that Assumption (ref) could be relaxed by requiring the bound $\mathcal{E}\left(\boldsymbol{\theta}\right)\geqslant c_{M}\Vert\boldsymbol{\theta}-\boldsymbol{\theta}_{0}\Vert_{2}^{2}$ to hold only for certain sparse vectors $\boldsymbol{\theta}$. We have opted for a less general statement to avoid additional technicalities.

The following assumption requires additional technical regularity of the loss function.

Assumption (ref).(ref) states that the loss function is locally Lipschitz in the first argument with the Lipschitz “constant” $L(\boldsymbol{x},\boldsymbol{y})$ being sufficiently well-behaved. The local Lipschitzness required in (ref) actually follows from the loss convexity in Assumption (ref) rockafellar_convex_1970, so Assumption (ref).(ref) should be regarded as a mild moment condition. Assumptions (ref).(ref) and (ref).(ref) essentially state that, viewed as functions of $\boldsymbol{\theta}$, both the loss and its derivative are mean-square continuous at $\boldsymbol{\theta}_0$. When the loss $m(\cdot,\boldsymbol{y})$ is globally Lipschitz uniformly in $\boldsymbol{y}$ (thus allowing the choice $c_L=\infty$), Assumption (ref).(ref) boils down to the regressors having sufficiently many absolute moments, and Assumption (ref).(ref) reduces to the requirement that the largest eigenvalue of $\mathrm{E}[\boldsymbol{X}\boldsymbol{X}^\top]$ is bounded from above.\footnote{Boundedness of eigenvalues is a standard assumption in the semi- and non-parametric estimation literature. See e.g. belloni_new_2015 and soerensen_2024.} Examples of globally Lipschitz losses are the logit likelihood loss in Example (ref) and the trimmed LAD loss in Example (ref).

This assumption is a sparsity condition, stating that $\boldsymbol{\theta}_{0}$ lies in an $\ell_{q}$-“ball” of “radius” $s_{q}^{1/q}$. We interpret the $q=0$ case in the limiting sense $\lim_{q\to0_{+}}\sum_{j=1}^{p}|\theta_{0,j}|^{q}=\sum_{j=1}^{p}\mathbf{1}(\theta_{0,j}\neq0)$ so as to nest the case of exact sparsity with (at most) $s_{0}$ non-zero entries. When $q>0$, we have only approximate sparsity, allowing possibly many---but typically small---non-zero entries. Related notions of sparsity appear in many papers on estimation of high-dimensional models. See Remark (ref) for further discussion.

Under Assumption (ref), we can (almost surely) define $\boldsymbol{S}_n\in\mathbb{R}^p$ by

In this paper we refer to $\boldsymbol{S}_n$ as the score.

We are now ready to present a theorem that provides probabilistic guarantees for $\ell_1$ and $\ell_2$ estimation errors of the $\ell_{1}$-ME. The proof, given in Appendix (ref), builds on arguments of belloni_l1-penalized_2011. Related statements appear also in van_de_geer_high-dimensional_2008, bickel_simultaneous_2009, and negahban_unified_2012, among others. Although we could not find the exact same version of the theorem in the literature, we make no claims of originality for these bounds and include the theorem for expositional purposes and in order to motivate our method for choosing the penalty parameter $\lambda$.

To state the theorem, recall that we denote $\eta_n=\sqrt{\ln(pn)/n}$.

This theorem motivates our choice of the penalty parameter $\lambda$. Specifically, it demonstrates that we want a level of regularization sufficient to overrule the score $(\lambda\geqslant c_{0}\|\boldsymbol{S}_n\|_{\infty})$ with high probability, without making the penalty “too large” $(\lambda>\overline{\lambda}_n).$ An interested reader can also find an analogue of Theorem (ref) for the post-$\ell_1$-ME in Appendix (ref), but the general principle for choosing $\lambda$ remains the same.

Bootstrapping after Cross-Validation

We next provide a method for choosing the penalty parameter which is broadly available yet amenable to theoretical analysis. We split the section into two parts. In Section (ref), we discuss a generic bootstrap method that allows for choosing the penalty parameter $\lambda$ under availability of some generic estimators $\widehat{U}_i$ of $U_i:=m'_1(\boldsymbol{X}_i^\top\boldsymbol{\theta}_0,\boldsymbol{Y}_i),i\in[n]$. In Section (ref), we show how to obtain suitable estimators $\widehat U_i$ via cross-validation. By analogy with linear mean regression, we refer to $U:=m_1'(\boldsymbol{X}^\top\boldsymbol{\theta}_0,\boldsymbol{Y})$ as the residual.\footnote{The linear mean model $\mathrm{E}[Y|\boldsymbol{X}]=\boldsymbol{X}^\top\boldsymbol{\theta}_0$ and (half) square loss imply $U=\boldsymbol{X}^\top\boldsymbol{\theta}_0-Y$. The name “residual” stems from $U$ agreeing with the deviation $Y-\mathrm{E}[Y|\boldsymbol{X}]$ from the mean up to a sign.}

Bootstrapping the Penalty Level

Suppose for the moment that residuals $U_{i}=m_1'\left(\boldsymbol{X}_{i}^{\top}\boldsymbol{\theta}_{0},\boldsymbol{Y}_{i}\right)$ are observable. In this case, we can estimate the $(1-\alpha)$-quantile $q_n(1-\alpha)$ of $\|\boldsymbol{S}_n\|_{\infty} = \| \mathbb{E}_{n}[U_i\boldsymbol{X}_i] \|_{\infty}$ via the Gaussian multiplier bootstrap.\footnote{Recall that $\boldsymbol{S}_n$ is well-defined a.s. under Assumption (ref). We omit the qualifier throughout this section.} To this end, let $\{e_i\}_{i=1}^{n}$ be independent standard normal random variables that are independent of the data $\{(\boldsymbol{X}_i,\boldsymbol{Y}_i)\}_{i=1}^{n}$. Given that $\mathrm{E}[U\boldsymbol{X}]=\mathbf 0_p$ under mild regularity conditions, the Gaussian multiplier bootstrap estimates $q_n(1-\alpha)$ by \[ \widetilde{q}_n\left(1-\alpha\right):=\left(1-\alpha\right)\text{-quantile of }\max_{1\leqslant j\leqslant p}\left|\mathbb{E}_{n}\left[e_{i}U_{i}X_{i,j}\right]\right|\text{ given }\{(\boldsymbol{X}_{i},\boldsymbol{Y}_i)\}_{i=1}^{n}. \] Under certain regularity conditions, $\widetilde q_n(1-\alpha)$ delivers a good approximation to $q_n(1-\alpha)$, even if the dimension $p$ of the vectors $\boldsymbol{X}_i$ is much larger than the sample size $n$. To see why this is the case, let $\boldsymbol{Z}:=(Z_1,\dotsc,Z_p)^\top$ be a centered random vector in $\mathbb{R}^{p}$ and let $\{\boldsymbol{Z}_i\}_{i=1}^n$ be independent copies of $\boldsymbol{Z}$. As established in chernozhukov_gaussian_2013,chernozhukov_central_2017, the random vectors $\{\boldsymbol{Z}_i\}_{i=1}^n$ satisfy the following high-dimensional versions of the central limit and Gaussian multiplier bootstrap theorems: If for some constant $b\in (0,\infty)$ and a non-random sequence $\widetilde{B}_n$ in $[1,\infty)$, possibly growing to infinity, one has \[ \min_{1\leqslant j\leqslant p}\mathrm{E}[Z_{j}^{2}]\geqslant b,\quad\max_{k\in\{1,2\}}\max_{1\leqslant j\leqslant p}\mathrm{E}\left[|Z_{j}|^{2+k}\right]/\widetilde{B}_n^k\leqslant 1\quad\text{and}\quad\mathrm{E}\Big[\max_{1\leqslant j\leqslant p}Z_{j}^{4}\Big]\leqslant \widetilde{B}_n^4, \] then there is a constant $C_b\in(0,\infty)$, depending only on $b$, such that

and, with probability approaching one,

where $\mathcal{A}_{p}$ denotes the collection of all (hyper)rectangles in $\mathbb{R}^{p}$, and $\boldsymbol{\mathcal{N}}_n$ is a centered Gaussian random vector in $\mathbb{R}^p$ with covariance matrix $\mathrm{E}[\boldsymbol{Z}\boldsymbol{Z}^{\top}]$. Provided $\widetilde{B}_n^4\ln^{7}(pn)/n\to0$, applying these two results with $\boldsymbol{Z}_i = U_i\boldsymbol{X}_i$ for all $i\in[n]$ and noting that sets of the form $A_t = \{\boldsymbol{u}\in\mathbb{R}^p;\max_{1\leqslant j\leqslant p}|u_j|\leqslant t\}$, $t\in[0,\infty)$, are indeed rectangles, suggest that the Gaussian multiplier bootstrap estimator $\widetilde q_n(1-\alpha)$ provides a good approximation to $q_n(1-\alpha)$.

As we typically do not observe the residuals $U_i = m_{1}'(\boldsymbol{X}_{i}^{\top}\boldsymbol{\theta}_{0},\boldsymbol{Y}_{i})$, the method described above is infeasible. Fortunately, the result ((ref)) continues to hold upon replacing $\{\boldsymbol{Z}_{i}\}_{i=1}^n$ with estimators $\{\widehat{\boldsymbol{Z}}_{i}\}_{i=1}^{n}$, provided these estimators are “sufficiently good,” in the sense to be defined below; see (ref). Suppose therefore that residual estimators $\{\widehat{U}_{i}\}_{i=1}^{n}$ are available. We then compute

and a penalty level follows as

We refer to this method for obtaining a penalty level as the bootstrap method (BM) and to $\widehat{\lambda}^{\mathtt{bm}}_{\alpha}$ itself as the bootstrap penalty level.

To ensure that $\widehat{q}^{\texttt{bm}}\left(1-\alpha\right)$ indeed delivers a good approximation to $q_n(1-\alpha)$, we invoke the following assumption.

This assumption imposes a few minor regularity conditions. It requires, in particular, that all components of the vector $\boldsymbol{X}$ are normalized to be on the same scale. Since this assumption is high level, we verify it under low-level conditions in Appendix (ref) for the examples in Section (ref).

Our next result provides convergence rates for the $\ell_1$-ME based on the bootstrap method.

The idea of using a bootstrap procedure to select the penalty level in high-dimensional estimation is in itself not new. chernozhukov_gaussian_2013 use a Gaussian multiplier bootstrap to tune the Dantzig selector candes_dantzig_2007 for the high-dimensional linear model allowing both non-Gaussian and heteroskedastic errors. Note, however, that chernozhukov_gaussian_2013 presumes access to a preliminary Dantzig selector, which is used to estimate residuals. The condition ((ref)) is similarly high-level in the sense that it does not specify how one performs residual estimation in practice. Our primary contribution lies in providing a method for coming up with good residual estimators. We turn to this task in the next subsection, where we also compare the rates with those appearing in the literature and discuss the side conditions under which they are derived; see Remarks (ref) and (ref).

Cross-Validating Residuals

In this subsection, we explain how residual estimation can be performed via cross-validation (CV). To describe our CV residual estimator, fix any integer $K\geqslant2$, and let $\{I_{k}\}_{k=1}^K$ partition the observation indices $[n]$. Provided $n$ is divisible by $K$, the even partition

is natural, but not necessary. For the formal results below, we only require that each $I_{k}$ specifies a “substantial” subsample; see Assumption (ref) below.

Let $\Lambda_{n}$ denote a finite subset of $(0,\infty)$ composed by candidate penalty levels. In Assumption (ref) below, we require $\Lambda_{n}$ to be “sufficiently rich.” Our CV procedure then goes as follows. First, estimate the vector of parameters $\boldsymbol{\theta}_{0}$ by

for each candidate penalty level $\lambda\in\Lambda_{n}$ and holding out each subsample $k\in[K]$ in turn. Second, determine the penalty level

by minimizing the out-of-sample loss over the set of candidate penalties. Third, estimate residuals $U_{i}=m_{1}'(\boldsymbol{X}_{i}^{\top}\boldsymbol{\theta}_{0},\boldsymbol{Y}_{i}),i\in[n]$, by predicting out of each estimation subsample, i.e.,

Note here that since $I_k$ and $I_k^c$ have no elements in common, the derivative $m_1'(\boldsymbol{X}^{\top}_i\widehat\boldsymbol{\theta}_{I_k^c}(\lambda),\boldsymbol{Y}_i)$ exists for all $i\in I_k$, $k\in[K]$, and $\lambda\in\Lambda_n$ almost surely by Assumption (ref). The residual estimates $\{\widehat U_i^{\mathtt{cv}}\}_{i=1}^n$ are therefore almost-surely well-defined even though the function $m(\cdot,\boldsymbol{y})$ is not necessarily differentiable.

Combining the bootstrap penalty level $\widehat{\lambda}^{\mathtt{bm}}_\alpha=c_{0}\widehat{q}^{\texttt{bm}}\left(1-\alpha\right)$ from the previous subsection with the CV residual estimates $\widehat{U}_i=\widehat{U}_{i}^{\mathtt{cv}}$ from this subsection, we obtain the {\em bootstrap-after-cross-validation} (BCV) method for estimating the quantile $q_n(1-\alpha)$,

and the BCV penalty level follows as

To analyze the $\ell_1$-ME implied by BCV, we invoke the following two assumptions.

Assumption (ref) means that we rely upon classical $K$-fold cross-validation with fixed $K$. This assumption does rule out leave-one-out cross-validation, since $K=n$ and $I_{k}=\left\{ k\right\} $ imply $|I_k|/n\to 0$. Assumption (ref) allows for a rather large candidate set $\Lambda_{n}$ of penalty values. Note that the largest penalty value, $C_{\Lambda}$, can be set arbitrarily large and the smallest value, $c_{\Lambda}/n$, converges rapidly to zero. As a part of the proof of Theorem (ref) below, we show that these properties ensure that the set $\Lambda_{n}$ eventually contains a “good” penalty candidate, say $\lambda_{\ast}$, in the sense of leading to a uniform bound on the excess risk of subsample estimators $\widehat{\boldsymbol{\theta}}_{I_{k}^{c}}\left(\lambda_{\ast}\right),k\in[K]$ and, because of that, the CV residual estimators are reasonable inputs for the bootstrap method, in the sense of satisfying (ref). Combining this finding with Lemma (ref), we obtain convergence rates for the $\ell_1$-ME implied by BCV.

We next consider the post-$\ell_1$-penalized M-estimator (post-$\ell_1$-ME). The main motivation for the post-$\ell_1$-ME is that the $\ell_1$-ME may be severely biased because it shrinks coefficients toward zero. By refitting the non-zero coefficients of the $\ell_1$-ME without the penalty in the criterion function in (ref), the post-$\ell_1$-ME attempts to reduce this bias.

To define the post-$\ell_1$-ME, for any $\overline\boldsymbol{\theta}\in\Theta$, we define the set $\widetilde\Theta(\mathrm{supp}(\overline\boldsymbol{\theta}))\subseteq\Theta$ by

Then for any $\ell_1$-ME, i.e. a solution $\widehat\boldsymbol{\theta}\in\widehat\Theta(\lambda)$ to the optimization problem in (ref), the corresponding post-$\ell_1$-ME is defined as any element $\widetilde\boldsymbol{\theta}$ of the set $\widetilde\Theta(\mathrm{supp}(\widehat\boldsymbol{\theta}))$. Note that there could be multiple post-$\ell_1$-MEs.\footnote{As in our treatment of $\ell_1$-ME, we implicitly assume that a post-$\ell_1$-ME exists.} Our treatment below covers the set $\widetilde\Theta(\lambda)$ of all post-$\ell_1$-MEs, which we denote

To analyze the post-$\ell_1$-ME, we will use the following two additional assumptions.

Assumption (ref) strengthens the almost-sure differentiability in Assumption (ref). The stronger smoothness requirement precludes the trimmed LAD loss function in Example (ref), but it can be easily verified under more low-level conditions for the trimmed LS loss function in the same example as well as for all other examples from Section (ref); see Appendix (ref). Assumption (ref) is satisfied if the entries of $\boldsymbol{X}$ are independent standard Gaussian, for example. Related assumptions appear in the existing literature on high-dimensional estimation.

With these added assumptions, we can derive the convergence rates for the post-$\ell_1$-ME.

Debiased Estimation and Inference

In this section, we describe how to construct $\sqrt n$-consistent and asymptotically normal estimators of individual components of the vector $\boldsymbol{\theta}_0$ defined in (ref). Since these estimators are asymptotically unbiased and have easily estimable asymptotic variance, they lead to standard inference procedures for testing hypotheses about and building confidence intervals for individuals components of $\boldsymbol{\theta}_0$. Our approach here is based on the concept of Neyman orthogonal equations and closely follows the literature on double/debiased machine learning chernozhukov2018double. We note that the tools developed in this section rule out the trimmed LAD loss in Example (ref), as this function is not sufficiently smooth to satisfy our Assumption (ref).

Without loss of generality, suppose that we are interested in the first component of the vector $\boldsymbol{\theta}_0$, so that $\boldsymbol{\theta}_0 = (\beta_0,\boldsymbol{\gamma}_0^{\top})^{\top}$, where $\beta_0\in\mathbb{R}$ is the scalar parameter of interest and $\boldsymbol{\gamma}_0\in\mathbb{R}^{p-1}$ is a vector of nuisance parameters. To derive a $\sqrt n$-consistent and asymptotically normal estimator of $\beta_0$, write $\boldsymbol{X} = (D,\boldsymbol{W}^{\top})^{\top}$, so that $\boldsymbol{X}^{\top}\boldsymbol{\theta}_0 = D\beta_0 + \boldsymbol{W}^{\top}\boldsymbol{\gamma}_0$, and let $\boldsymbol{\mu}_0\in\mathbb{R}^{p-1}$ be a vector that is defined as a solution to the following system of equations:

Note that this system has a solution $\boldsymbol{\mu}_0$ and this solution is unique as long as the matrix $\mathrm{E}[m''_{11}(\boldsymbol{X}^{\top}\theta_0,\boldsymbol{Y})\boldsymbol{W}\boldsymbol{W}^{\top}]$ is non-singular, which is the case under our assumptions.\footnote{See Lemma (ref) in the appendix for the precise statement.} With this definition in mind, by inspecting the first-order conditions associated with (ref), we have

We obtain an estimator of $\beta_0$ by solving an empirical version of this equation, where we replace the (high-dimensional) vectors $\boldsymbol{\gamma}_0$ and $\boldsymbol{\mu}_0$ by suitable estimators. Here, $\sqrt n$-consistent and asymptotically normal estimation of $\beta_0$ is possible due to (ref) being Neyman orthogonal with respect to $\boldsymbol{\gamma}_0$ and $\boldsymbol{\mu}_0$, which means that this equation is first-order insensitive with respect to perturbations in $\boldsymbol{\gamma}_0$ and $\boldsymbol{\mu}_0$. Specifically, we have

which follows from (ref) and (ref), respectively. Neyman orthogonality thus facilitates simple inference for the low-dimensional $\beta_0$ despite possibly complicated estimation of the high-dimensional $\boldsymbol{\gamma}_0$ and $\boldsymbol{\mu}_0$. Formally, we consider the following procedure:

Note that (even without refitting) this procedure gives two estimators of $\beta_0$: $\widetilde\beta$ on the first step and $\widehat\beta$ on the third step. As it turns out, the estimator $\widehat\beta$ is better, in the sense that it can be established as both asymptotically unbiased, $\sqrt n$-consistent and asymptotically normal. To derive these properties, we impose the following assumptions.

Assumption (ref) essentially means that there is non-trivial variation in the variable of interest $D$ after partialling out the controls $\boldsymbol{W}$. In the familiar case of the linear mean model with square loss, non-trivial variation follows from the usual rank condition for identification of $\boldsymbol{\theta}_0$; see Appendix (ref) for details.\footnote{More generally, Assumption (ref) is implied by the eigenvalues of the matrix $\mathrm{E}[U^2\boldsymbol{X} \boldsymbol{X}^{\top}]$ being bounded away from zero, which is a non-degeneracy condition.} Assumption (ref) imposes minor regularity conditions requiring a certain amount of integrability of the random variables in the model and transformations thereof.

For inference purposes, we also invoke a stronger smoothness condition.

This assumption strengthens Assumption (ref) from Section (ref) (which is why we reuse the symbol $C_m$ for the constant). Note that Assumption (ref) does not hold for the trimmed LAD loss function in Example (ref), which means that our inference approach does not apply for this loss function. In addition, Assumption (ref) does not hold for the trimmed LS loss function in the same example whenever $\boldsymbol{\theta}_0 = \mathbf 0_{p}$. Although we believe it should be possible to perform inference in these cases using methods from BCK17 developed for the case of a high-dimensional linear quantile regression model, we leave this line of work for the future. In Appendix (ref), we verify Assumption (ref) for all other examples from Section (ref) including Example (ref) with the trimmed LS loss function whenever $\boldsymbol{\theta}_0 \neq \mathbf 0_p$.

The next assumption controls the impact of points of non-smoothness in the loss, if any.

Assumption (ref) holds trivially in Examples (ref) and (ref), as for those examples Assumption (ref) holds with $J=1$. When combined with the requirement that $\overline\Delta_n\to0$ (sufficiently fast), Assumption (ref) does impose quite a bit of structure in Examples (ref) and (ref), however. In Example (ref), this assumption is satisfied if the conditional distribution of $Y$ given $\boldsymbol{X}$ is absolutely continuous with bounded PDF. In Example (ref) with the trimmed LS loss function, it is satisfied if the conditional distribution of $Y_1$ given $(\boldsymbol{X},Y_1>0)$, the conditional distribution of $Y_2$ given $(\boldsymbol{X},Y_2>0)$, and the (unconditional) distribution of $\boldsymbol{X}^{\top}\boldsymbol{\theta}_0$ are all absolutely continuous with bounded (uniformly over $n$) PDFs; see Appendix (ref) for details.\footnote{Note here that the requirement that the distribution of $\boldsymbol{X}^{\top}\boldsymbol{\theta}_0$ is absolutely continuous with bounded PDF implies that $\boldsymbol{\theta}_0$ is sufficiently well separated from $\mathbf 0_p$.}

Assumption (ref) is a high-level assumption placed on the estimators from Steps 1 and 2 in Algorithm (ref). When $\widetilde\boldsymbol{\theta}$ is $\ell_1$-ME or post-$\ell_1$-ME based on BCV, we can lean on the bounds from Theorem (ref). For the estimation error of $\widetilde\boldsymbol{\mu}$, however, we cannot use Theorem (ref), as this estimator does not fit into our framework because of the presence of estimated weights in the optimization problems (ref) and (ref). However, these optimization problems correspond to LASSO and post-LASSO with estimated weights, and such estimators are well studied in the literature. See e.g. belloni2016post, where one can find the appropriate rates for the estimation error of $\widetilde\boldsymbol{\mu}$ in terms of the sparsity of $\boldsymbol{\mu}_0$.

We next present a theorem on the asymptotic distribution of the debiased estimator $\widehat\beta$.

This theorem shows that the estimator $\widehat\beta$ is asymptotically unbiased and normal under plausible regularity conditions. The “asymptotic” variance $\sigma_0^2$ appearing in this theorem, which depends on $n$ in general via the distribution $P$ of $(\boldsymbol{X},\boldsymbol{Y})$, is easily estimable. For example, one can use a plug-in estimator

with the estimators $\widetilde\beta,\widetilde\boldsymbol{\gamma}$ and $\widetilde\boldsymbol{\mu}$ stemming from Steps 1 and 2 of Algorithm (ref) (possibly with refitting) using BCV as the penalty rule in both steps. Alternatively, one can incorporate Step 3 of the same algorithm and use

It is rather standard to derive consistency of these estimators. Also, because of asymptotic normality of $\widehat\beta$, it is then straightforward to perform inference on $\beta_0$. For example, an asymptotically valid $(1-\alpha)\times100\%$ confidence interval for $\beta_0$ takes the standard form $[\widehat\beta - z_{\alpha/2}\widehat\sigma/\sqrt n,\widehat\beta + z_{\alpha/2}\widehat\sigma/\sqrt n]$, where $\widehat\sigma$ is given either by (ref) or by (ref), and $z_{\alpha/2}$ is the $(1-\alpha/2)$-quantile of the standard normal distribution.

Simulations

In this section we investigate the finite-sample behavior of our estimators based on the bootstrap-after-cross-validation (BCV) method for obtaining penalty levels proposed in Section (ref). We also compare our estimation and inference methods to ($K$-fold) cross-validation, which lacks general theoretical justification but is a popular method in practice.

Simulation Design

We consider a master data-generating process (DGP) of the form \[ Y_{i}=\mathbf{1}\left(\beta_{0}D_{i}+\sum_{j=1}^{p-1}\gamma_{0j}W_{i,j}+\varepsilon_{i}>0\right),\quad\varepsilon_{i}\mid D_{i},\boldsymbol{W}_{i}\sim\mathrm{N}\left(0,1\right),\quad i\in\left[n\right], \] thus implying a binary probit model as in Example (ref). The regressors $\boldsymbol{X}=(D,\boldsymbol{W}^{\top})^{\top}$ are distributed jointly centered Gaussian $\boldsymbol{X}\sim\mathrm{N}(\boldsymbol{0},\boldsymbol{\Sigma}(\rho))$ with covariances (and correlations) \[ \Sigma_{j,k}(\rho):=\mathrm{cov}\left(X_j,X_k\right)=\mathrm{E}\left[X_jX_k\right]=\rho^{\left|j-k\right|},\quad\left(j,k\right)\in\left[p\right]^2. \] Hence, the regressor covariance matrix $\boldsymbol{\Sigma}\left(\rho\right)$ takes a Toeplitz form with the overall correlation level being dictated by $\rho$. We allow $\rho\in\left\{ 0,.2,\dotsc,.8\right\} $, thus running the gamut of (positive) correlation levels. Since $\varepsilon_{i}$'s are standard normal, the “noise” $\mathrm{var}\left(\varepsilon\right)$ in our DGP is fixed at one. Hence, the signal-to-noise ratio (SNR) equals the “signal,” \[ \mathrm{SNR}:=\frac{\mathrm{var}(\boldsymbol{X}^{\top}\theta_{0})}{\mathrm{var}(\varepsilon)} =\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}\left(\rho\right)\boldsymbol{\theta}_{0}. \] which depends on both the correlation level and coefficient pattern. We consider the patterns:

The exactly sparse pattern has only non-zero coefficients for the first couple of regressors $(s_{0}=2)$, and both non-zeros are clearly separated from zero, thus allowing perfect variable selection. The implied signals (hence SNRs) are

Compared to existing simulation studies for high-dimensional binary response models, the SNRs considered here are relatively low.\footnote{For example, the binary logit designs in friedman_regularization_2010 and ng_feature_2004 imply SNRs of three and over 30, respectively.}

Note that the SNR is increasing with the regressor correlation, such that sampling from a high-$\rho$ DGP tends to produce an easier estimation problem compared to sampling from a low-$\rho$ DGP, keeping all other things equal. When reporting results below for $\rho=0$ (our baseline), we are thus considering the worst correlation scenario.\footnote{The same comments apply to the other coefficient patterns albeit with the more complicated signal \[ \mathrm{var}\left(\boldsymbol{X}^{\top}\boldsymbol{\theta}_{0}\right)=\sum_{j=1}^{p}\theta_{0,j}^{2} +2\sum_{j=1}^{p-1}\sum_{k=j+1}^{p}\theta_{0,j}\theta_{0,k}\rho^{k-j}. \] }

In contrast to the exactly sparse pattern, the approximately sparse pattern involves all non-zeros $(s_{0}=p)$, which are not bounded away from zero, such that variable selection mistakes are bound to happen. To see that this pattern is in fact approximately sparse, note that for every $q\in(0,1]$ one has $ \sum_{j=1}^{p}\left|\theta_{0,j}\right|^{q}\leqslant\sum_{j=1}^{\infty}\left|\theta_{0,j}\right|^{q}=1/(1-2^{-q/2}). $ Hence, for the purpose of Assumption (ref), we can choose $q\in(0,1]$ freely and pair it with $s_q=1/(1-2^{-q/2})$. The base $1/\sqrt{2}$ of the approximately sparse pattern was here chosen to (approximately) equate the signals arising from the approximately and exactly sparse coefficient patterns in the baseline case of uncorrelated regressors $(\rho=0)$, which amounts to $\|\boldsymbol{\theta}_{0}\|_{2}^{2}$. The relevance of a regressor, as measured by its coefficient, is rapidly decaying in the regressor index $j$, such that the vast majority of the signal is captured by a small fraction of the regressors. For example, in the baseline case of uncorrelated regressors $(\rho=0)$, the first 10 regressors account for 99.9 percent of the signal (two).

In between these two extremes lies the intermediate pattern. This pattern was created by cutting off the approximately sparse coefficient sequence at the smallest regressor index $j^{\ast}$, such that regressors $[j^{\ast}]$ account for at least 95 percent of the baseline signal. (Here: $j^{\ast}=5$.) For this pattern, perfect variable selection is possible but unlikely.

We consider sample sizes $n\in\left\{ 100,200,400\right\} $ and limit attention to the high-dimensional regime by fixing $p=n$ throughout.

Estimation and Implementation

We consider the following four estimators arising from $\ell_1$-ME (ref) and post-$\ell_1$-ME (ref) based on either the CV or BCV penalty levels in (ref) and (ref), respectively:

When discussing normal approximations based on three-step debiasing (Algorithm (ref)), we use the same method in both Steps 1 and 2. For example, the “post-BCV” inference procedure refers to post-BCV in the first step, followed by post-BCV in the second step (i.e. both optional steps are taken).

Our BCV and post-BCV estimation methods require us to specify a score markup $c_{0}\in(1,\infty)$ and probability tolerance rule $\alpha=\alpha_{n}$. We here follow the recommendation in belloni_sparse_2012 for the LASSO and post-LASSO and take $c_{0}=1.1$ and $\alpha_{n}=.1/\ln(p\lor n)$ as our benchmark. The latter function, slowly decaying in $p\lor n$, leads to $\alpha\approx2.2\%,1.9\%$ and $1.7\%$ for $n=100,200$ and $400$, respectively. We also look at the alternative score markups $\{1,1.05\} $, the first one being excluded by the theory in Section (ref). The alternative probability tolerance rule $\alpha_{n}=10/n$ leads to qualitatively identical conclusions, cf. Appendix (ref). We stress that the benchmark choices of $c_0$ and $\alpha_n$ are only rules of thumb that tend to perform well in the simulation designs considered here. Other choices of score markups and probability tolerance rules may have better properties in other DGPs.

We have previously treated all coefficients in the same manner, in that they are all penalized and with equal weight. However, in an empirical application one is typically confident that an intercept belongs in the model. For this reason, the (intercept) coefficient on the constant regressor is usually not penalized during estimation. Moreover, to justify equal penalty weighting, prior to estimation one typically brings the (non-constant) regressors onto the same scale by dividing them by their respective sample standard deviations. To align our simulation study with these empirical practices, we include unpenalized intercepts in both Steps 1 and 2 of Algorithm (ref) and rescale regressors. (The intercepts are still suppressed in our notation.) That is, we treat neither the zero (true) intercept nor equivariant regressors as information known to the researcher. In these aspects our simulations are therefore empirically calibrated.

For each sample size $n(=p)$, each correlation level $\rho$, and each coefficient pattern, we use 2,000 independent simulation draws and 1,000 independent standard Gaussian bootstrap draws per simulation draw and per estimation step (when relevant). We assign observations to $K$ approximately equally large folds $\{ I_{k}\} _{k=1}^{K}$ for both the first and second steps, shuffling the assignments in between. We keep $K=3$ throughout and use the same folds for all estimators to facilitate comparison.\footnote{Three folds is the minimum value allowed by cv.glmnet. Preliminary and unreported simulation experiments suggest that using 5-fold (instead of 3-fold) CV only affects the average errors reported below at the third decimal. Similarly, using 2,000 Gaussian bootstraps (instead of 1,000) appears to only affect these averages at the fourth decimal.}

All simulations are carried out in R with cross-validation done using glmnet::cv.glmnet, and refitting done using stats::glm.\footnote{We use R version 4.2.2 and glmnet version 4.1-6.} When constructing the candidate penalty set $\Lambda_{n}$, we use the glmnet default settings, which creates a log-scale equi-distant grid of a 100 candidate penalties from the threshold penalty level to essentially zero. The threshold is the (approximately) smallest level of penalization needed to set every coefficient to zero, thus resulting in a trivial (null) model.\footnote{Log-scale equi-distance from a “large” candidate value to essentially zero fits well with the form of $\Lambda_{n}$ in our Assumption (ref) (interpreting $c_{\Lambda}/n\approx0$). However, the threshold penalty is a function of the data and, thus, random. The resulting candidate penalty set used in our simulations is therefore also random, and thus, strictly speaking, not allowed by Assumption (ref). Moreover, the number of candidate values $|\Lambda_{n}|$ is here held fixed. We believe these deviations from our theory to be only a minor issue.}

Note that cv.glmnet calculates and stores the out-of-fold linear forms $\boldsymbol{X}_{i}^{\top}\widehat{\boldsymbol{\theta}}_{I_{k}^{c}}(\lambda)$ (with an intercept, if relevant) for each $i\in I_{k}$, fold $k$ and candidate penalty $\lambda$, and allows for extraction of estimates for penalty levels off the regularization grid via linear interpolation. Hence, compared to CV, there is essentially zero added computational burden associated with using BCV.

Results

Non-Existence and Treatment of Missing Values

While the $\ell_{1}$-penalized probit estimators BCV and CV always exist (cf. Section (ref)), refitting after variable selection based on either of these estimator can fail. For example, in our binary response setting, without any penalty one may encounter complete separation of the outcomes based on the fitted probabilities, in which case the refitted estimates fail to exist (as real numbers).\footnote{Strictly speaking, the $\ell_{1}$-penalized probit estimator fails to exist when all outcomes are of the same label and some coefficient (here: the intercept) goes unpenalized. In none of our simulated datasets did we encounter all zeros or all ones. See Appendix (ref) and, in particular, Remark (ref) for more discussion.}

Across all simulation designs and draws, refitting after CV fails in nearly 15% of all cases. The fraction of such non-existent post-CV cases varies with the DGP and can be higher than 47%. Since post-CV estimation and debiasing procedures do not appear well-defined in our context, we drop them from further consideration.

In contrast, refitting after BCV fails to converge in only about 0.01% of all cases.\footnote{Specifically, convergence fails in 74 out of a total of 540,000 cases, where the total equals the product of the numbers of simulation draws (2,000), correlation levels (5), sample/problem sizes (3), coefficient patterns (3), score markups (3) and probability tolerance rules (2).} Since we find this fraction miniscule, when reporting results below we choose to simply omit the problematic cases from the relevant post-BCV statistics; see also the figure notes.

Estimation Error

Figure (ref) shows the mean $\ell_{2}$ estimation error (for the slope coefficients, averaging over the 2,000 simulation draws) arising from BCV, post-BCV and CV, respectively, using benchmark tuning.

Mean $\ell_{2}$ estimation error is here depicted as a function of the sample/problem size (the tile column), coefficient pattern (the tile row), and correlation level (the horizontal axis in each tile). The horizontal line at $\|\boldsymbol{\theta}_{0}\|_{2}$ facilitates comparison with the trivial “estimator” $\widehat{\boldsymbol{\theta}}\equiv\boldsymbol{0}_p$.

One observation evident from this figure is that the error curves of BCV, post-BCV and CV can cross. Hence, these estimators cannot be ranked in terms of mean $\ell_{2}$ estimation error, in general. However, for the largest sample/problem size considered, post-BCV outperforms CV for small to medium levels of correlation, and CV outperforms BCV.\footnote{We reach qualitatively identical conclusions from inspecting the median $\ell_{2}$ estimation errors. Hence, these findings are not limited to one particular feature of the error distributions. We omit the corresponding median plots due to their similarity with the mean error plots. Figures are available upon request.}

Increasing the sample size (moving from left to right) leads to a downward shift in mean estimation error for all three estimators, which is indicative of convergence. Convergence appears to take place no matter the coefficient pattern or regressor correlation level even though the number of candidate regressors matches the sample size. Increasing the number of non-zeros in $\boldsymbol{\theta}_{0}$ (moving from top to bottom) leads to an upward shift in mean estimation error. This finding is consistent with convergence slowing down with $q$ and $s_q$ as predicted by Theorems (ref) and (ref).

We next investigate the impact of the choice of score markup $c_{0}$. Figures (ref) and (ref) plot the mean $\ell_{2}$ estimation error for $c_{0}=1,1.05$ and (the previously used) $1.1$, each sample/problem size and coefficient pattern, and for the BCV and post-BCV estimators, respectively.

Figure (ref) suggests that increasing $c_{0}$ away from one slightly worsens (mean estimation error) performance of BCV. While our theory takes $c_{0}$ strictly greater than one, any value near one---including the limit case of one itself--- appears to lead to near identical results.\footnote{That mean BCV error is downward sloping for small to moderate $\rho$ levels is due to the signal being increasing in $\rho$ and need not translate to other correlation or coefficient patterns.} For post-BCV (Figure (ref)), the findings are similar. In fact, at least for the largest sample/problem size, the exact value of $c_{0}\in\left\{ 1,1.05,1.1\right\} $ has little to no impact on mean error. Note that our findings for post-BCV apply even with the approximately sparse coefficient pattern, where variable selection mistakes are bound to occur.

To conclude this subsection, we note that it is a well-known puzzle in the LASSO literature that the theory typically requires that $c_0$ is strictly bigger than one, with the estimation error bounds deteriorating as $c_0$ approaches one, while simulation experience suggests that the estimation errors are insensitive with respect to $c_0$ when $c_0$ is close to one. We believe that solving this puzzle remains one of the key challenges in this literature.

Normal Approximation

We next assess the normal approximations resulting from three-step debiasing (Algorithm (ref)) using either BCV, post-BCV or CV. Instead of looking at the standardized estimate $\sqrt{n}(\widehat{\beta}-\beta_{0})/\sigma_{0}$ for the true asymptotic variance $\sigma_0^2$ given in (ref), we form an estimate $\widehat{\sigma}^{2}$ and consider the studentized estimate $\sqrt{n}(\widehat{\beta}-\beta_{0})/\widehat{\sigma}$. That is, we take into account the unknown nature of the $\sigma_{0}^{2}$, as required in an empirical application.

To construct the estimate $\widehat{\sigma}^{2}$, we first leverage the binary response model to establish the (conditional information) equality $ \mathrm{E}[m_{1}'\left(\boldsymbol{X}^{\top}\boldsymbol{\theta}_{0},\boldsymbol{Y}\right)^{2}\mid \boldsymbol{X}]=\mathrm{E}[m_{11}''\left(\boldsymbol{X}^{\top}\boldsymbol{\theta}_{0},\boldsymbol{Y}\right)\mid \boldsymbol{X}]. $ We then use the definition of $\boldsymbol{\mu}_{0}$ to establish the (weighted projection) equality \[ \mathrm{E}\big[m_{11}''\left(\boldsymbol{X}^{\top}\boldsymbol{\theta}_{0},\boldsymbol{Y}\right)\left(D-\boldsymbol{W}^{\top}\boldsymbol{\mu}_{0}\right)^2\big]=\mathrm{E}\left[m_{11}''\left(\boldsymbol{X}^{\top}\boldsymbol{\theta}_{0},\boldsymbol{Y}\right)\left(D-\boldsymbol{W}^{\top}\boldsymbol{\mu}_{0}\right)D\right]. \] Again using the binary response model, we can evaluate

where $f$ and $F$ denote the PDF and CDF, respectively, associated with the binary response model.\footnote{In our current binary probit setting, these functions are the standard normal PDF and CDF, respectively. In the empirical application in Section (ref), we also use the logistic distribution, leading to the binary logit.} This allows us to simplify the expression for $\sigma_{0}^{2}$ to \[ \sigma_{0}^{2}=1\big\slash\mathrm{E}\left[\omega_{F}\left(\boldsymbol{X}^{\top}\boldsymbol{\theta}_{0}\right)\left(D-\boldsymbol{W}^{\top}\boldsymbol{\mu}_{0}\right)D\right] \] and leads to the following estimator of $\sigma_0^2$: \[ \widehat{\sigma}^{2}:=1\big\slash\mathbb{E}_{n}\big[\omega_{F}\big(\widehat{\beta}D_{i}+\boldsymbol{W}_{i}^{\top}\widetilde{\boldsymbol{\gamma}}\big)\big(D_{i}-\boldsymbol{W}_{i}^{\top}\widetilde{\boldsymbol{\mu}}\big)D_{i}\big], \] with $\widetilde{\boldsymbol{\gamma}},\widetilde{\boldsymbol{\mu}}$ and $\widehat{\beta}$ given by Steps 1, 2 and 3, respectively, of Algorithm (ref) given different rules for choosing the penalties $\lambda_1$ and $\lambda_2$.\footnote{Alternatively, one can use the “sandwich” estimators (ref) and (ref). Experimenting with these estimators, we obtained numerically similar results as reported below for the estimator $\widehat{\sigma}^2$. We prefer $\widehat{\sigma}^2$ since it leverages both the binary response and projection structure.}

Figure (ref) shows the (kernel) densities of the studentized estimates using benchmark tuning and $\rho=0$.\footnote{All kernel densities are created using the R package ggplot2 with geom_density. In expectation of an approximately normal distribution, we use a Gaussian kernel and the silverman_density_1986 rule-of-thumb bandwidth (both geom_density defaults).}

The densities arising from BCV, post-BCV and CV, respectively, are here depicted as columns of tiles, where each tile row corresponds to a coefficient pattern and each graph within a tile a sample/problem size. Starting with the exactly sparse coefficient pattern (the top row), we see that both BCV and CV lead to considerable shrinkage bias even after debiasing the initial estimate of the focal parameter $\beta_{0}$. This feature is seen from the leftward shifts in the resulting densities compared to the standard normal density, here represented by the dotted line. These biases do not seem to disappear as $n$ increases, holding $n=p$. If anything, these distributions shift further left, which indicates that BCV requires a larger sample size. In constrast, the post-BCV density essentially collapses to the standard normal one, at least for $n=p=200$ and $400$.

As the coefficient pattern becomes less and less sparse (moving down), all approximations deteriorate, as is to be expected. While imperfect, the post-BCV densities are still decent approximations to the normal for both the intermediate and approximately sparse coefficient patterns. Moreover, only these densities appear to approach the standard normal as the sample/problem size increases.

While Figure (ref) depicts the normal approximations for the worst-correlation case $\rho=0$, in Figure (ref) we display the normal approximations as a function of $\rho$.

We here focus on the largest sample size $n(=p)=400$ and the (more challenging) approximately sparse coefficient pattern, again using benchmark tuning. Since post-BCV and CV appear to lead to better normal approximations than BCV, we display only results from the former two methods. Post-BCV leads to a relatively accurate normal approximation for every correlation level considered. Moreover, while the normal approximation stemming from CV appears to improve as $\rho$ increases, at no correlation level considered does CV lead to a visually better approximation than post-BCV.\footnote{We also investigated the robustness of the post-BCV-resulting normal approximations with respect to the markup $c_{0}\in\{1,1.05,1.1\}$. Parallelling our findings for mean estimation error in Figure (ref), the exact markup value appears to make little difference. Figures are available upon request.}

Revisiting Racial Differences in Police Use of Force

In this section we revisit the empirical setting in fryer_jr_empirical_2019 (henceforth: Fryer), who explored racial differences in police use of force. We here focus on the part of Fryer's regression analysis invoking the full Police-Public Contact Survey (PPCS) dataset with the outcome being an indicator for any use of force by the police (conditional on an encounter), thus leading to a binary response model as in our Example (ref).\footnote{See fryer_jr_empirical_2019 and the associated online appendix for alternative outcome variables and data sources as well as a detailed discussion of their relative merits and drawbacks.} Specifically, Fryer estimates models of the form

where $\mathrm{Force}$ indicates whether any force was used by the police when encountering a civilian, $\mathbf{Race}=(\mathrm{Black},\mathrm{Hisp},\mathrm{Other})^{\top}$ indicate the race of the civilian (black, hispanic and other than white, with white being the reference race), and $\boldsymbol{W}$ is a list of control variables (including a constant) capturing both civilian (e.g. gender and age), officer (e.g. majority race) and encounter characteristics (e.g. whether the civilian disobeyed, resisted or otherwise misbehaved).\footnote{See the fryer_jr_empirical_2019 notes for the full variable list and his online appendix A for descriptions.} Here $F$ is a placeholder for a strictly increasing known cumulative distribution function (CDF), which Fryer takes to be the logistic CDF $\Lambda$, thus leading to the binary logit model.

The PPCS logistic regression results reported in fryer_jr_empirical_2019 show that black and hispanic subjects are statistically significantly more likely to experience some form of force in interactions with the police, controlling for context and civilian behavior. We here look into the robustness of this finding by employing also

For brevity, we here single out the black dummy $(\mathrm{Black})$ and its coefficient $\beta_{\mathrm{Black}}$ and group the other non-white dummies ($\mathrm{Hisp}$ and $\mathrm{Other}$) with the controls (thus recasting $\boldsymbol{W}$ and $\boldsymbol{\gamma}_0$). We interpret the statement “there are no racial differences in police use of force” as there being no difference in the probability of force being used for black civilians relative to white subjects, holding everything else equal. That is, \[ \mathrm{P}\left(\mathrm{Force}=1\mid\mathrm{Black}=1,\boldsymbol{W}=\boldsymbol{w}\right) =\mathrm{P}\left(\mathrm{Force}=1\mid\mathrm{Black}=0,\boldsymbol{W}=\boldsymbol{w}\right) \] for all realizations $\boldsymbol{w}$ of the controls (with $\mathrm{Hisp}$ and $\mathrm{Other}$ both zero). Since the race dummies enter the strictly increasing $F$ in (ref) in an additive manner, using the model, such a zero probability difference is equivalent to a zero coefficient on the dummy for being black, i.e. $\beta_{\mathrm{Black}}=0$. We therefore take the latter as the hypothesis to be tested.

To this end, we first use the fryer_jr_empirical_2019 supplementary files and descriptions in his online appendix to recollect and recreate the PPCS dataset. Using the same supplementary files, we then replicate the PPCS logistic regression results in fryer_jr_empirical_2019 to all reported digits, which leaves us confident that we are indeed considering the original dataset.

We next apply three-step debiasing (Algorithm (ref)) with the loss function being either the negative logit or probit log-likelihood. Our simulation findings indicate that post-BCV debiasing outperforms both the BCV and CV equivalents. We therefore only consider the former.\footnote{For both Steps 1 and 2, we here use 10-fold cross-validation, i.e. $K=10$.} We use two sets of regressors. The first set (Basic Controls) corresponds to that in fryer_jr_empirical_2019, and is Fryer's largest set of controls. The only difference is that we include categorical regressors via dummies for their different levels, leaving one reference category for each. The second set of regressors (Basic Controls + Interactions) builds on the first by adding all first-order pairwise interactions between the controls (excluding the race dummies $\mathrm{Hisp}$ and $\mathrm{Other}$). After eliminating variables with zero variance or perfect correlation, the two sets include 30 and 327 non-constant regressors, respectively, which should be compared to a sample size of $n=59,668$ civilian--police encounters.\footnote{We note in passing that the 59,668 equals the total number of police-civilian encounters in the PPCS dataset covering the six surveys 1996, 1999, 2002, 2005, 2008 and 2011. The number of complete cases with respect to the regressors used is 9,930 and only leaves the years 2002 and 2011. For a clean comparison, we follow Fryer's approach to missing values.}

Table (ref) displays the $t$-values associated with testing the null hypothesis using either unpenalized or $\ell_{1}$-penalized methods.

Using only basic controls (the logit case being covered in Fryer), the $t$-statistics take on similar values for both unpenalized maximum likelihood (ML) and post-BCV methods. Thus, with only 30 non-constant candidate regressors, regularization has little impact. In contrast, the specification including both basic controls and interactions thereof leads to complete separation in the data, such that the (unpenalized) maximum likelihood estimates do not exist (as real numbers). For this case, some regularization is necessary. Hence, although the numbers of regressors considered here may not appear overwhelmingly large when compared to the sample size, the set of regressors is of great importance. Even when we include all first-order interactions between the controls, the $t$-statistics resulting from our three-step debiasing procedure remain of the same order as before.\footnote{The increase in the $t$-values for post-BCV upon inclusion of interactions is for both the logit and probit loss due to both a somewhat larger point estimate and a somewhat smaller standard error. The values underlying Table (ref) thus illustrate that more candidate regressors need not lead to larger standard error.} The $t$-tests based on our post-BCV debiasing lead us to reject the null hypothesis of no racial differences in police use of force at any reasonable significance level. This conclusion in fryer_jr_empirical_2019 therefore appears robust to the choice of controls.

To gauge the economic impact of our change in estimation procedures, we estimate the average partial effect (APE) of changing the civilian race from white to black. Iterating expectations and using ((ref)), the APE can be expressed as the average probability difference

where we bring back the other (non-white) civilian race dummies to clarify the comparison made. We estimate this APE by $ \widehat{\mathrm{APE}}_{\mathrm{Black}}:=\mathbb{E}_{n}[F\big(\widehat{\beta}_{\mathrm{Black}}+\boldsymbol{W}_{i}^{\top}\widehat{\boldsymbol{\gamma}}\big)-F\left(\boldsymbol{W}_{i}^{\top}\widehat{\boldsymbol{\gamma}}\right)], $ for point estimates $\widehat{\beta}_{\mathrm{Black}}$ and $\widehat{\boldsymbol{\gamma}}$ of $\beta_{\mathrm{Black}}$ and $\boldsymbol{\gamma}_{0}$, respectively. When results stem from (unpenalized) ML, we use the ML estimates. When results stem from three-step post-BCV debiasing, we use the debiased third-step estimate $\widehat{\beta}_{\mathrm{Black}}$ and the (biased) first-step estimate $\widehat{\boldsymbol{\gamma}}$. Table (ref) reports the APE estimates (in percentage points) corresponding to these procedures including either basic controls or basic controls with interactions.

Using only basic controls, (unpenalized) ML and post-BCV lead to APE estimates in the range of 1.1--1.4 percentage points regardless of the CDF used. (For context, the unconditional average of contacts in which PPCS respondents reported any force being used for white civilians is .7 percent.) Including interactions of basic controls, the post-BCV APE estimates roughly double in size to about 3 percentage points. Of course, these relatively large APE estimates may come with relatively large estimation error. However, as the APE in ((ref)) depends on many coefficients, it remains a non-trivial task to assign standard errors to these point estimates---a task falling outside the scope of this paper.

Finally, to get a feel for the computational burden associated with the methods proposed in this paper when applied to real data, in Table (ref) we report the computing time used by the above-mentioned estimation routines.

With only basic controls, the three-step post-BCV debiasing procedure takes at least ten times as long as (unpenalized) ML. This is not surprising, as the former method involves two rounds of (10-fold) CV, bootstrapping and refitting---no task of which is undertaken by ML. However, including also interactions, the ranking of the two approaches is reversed. The about tenfold increase in number of controls increases the computing time associated with post-BCV logit debiasing approximately linearly. For post-BCV probit debiasing, the corresponding increase is almost five fold. In contrast, as the ML estimates are not real numbers, without proper checks for solution existence, any (gradient-based) optimizer would iterate indefinitely in search of the ML estimates. We represent the non-existence of an ML estimate by infinite computing time.\footnote{While infinity may appear overly dramatic, we warn that glmnet does not check for optimizer divergence friedman_regularization_2010. We therefore opted for stats::glm for ML estimation and refitting.}

Of course, the Table (ref) runtimes are only single observations arising from our particular R implementation of our procedures, using a specific dataset, and our specific computing environment. As such, they need not translate to other settings.