The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
133,242 characters
Selecting Penalty Parameters of High-Dimensional M-Estimators using Bootstrapping after Cross-Validation
\title{Selecting Penalty Parameters of High-Dimensional M-Estimators using Bootstrapping after
Cross-Validation\thanks{Parts of this paper were previously circulated under the title ``Analytic
and Bootstrap-after-Cross-Validation Methods for Selecting Penalty Parameters of High-Dimensional
M-Estimators.'' We thank Richard Blundell, Victor Chernozhukov, Bo Honor\'{e}, Whitney Newey, Joris
Pinkse, Simon Reese, Azeem Shaikh, Mikkel S{\o}lvsten, Sara van de Geer and numerous seminar
participants for their insightful comments and discussions. Bohdan Salahub and Andrei Voronin
provided excellent research assistance. Chetverikov's work was supported by NSF Grant SES -
1628889.}} \author{Denis Chetverikov\footnote{Department of Economics, UCLA; e-mail:
\tt{[email removed]}.} \and Jesper R.-V. S{\o}rensen\footnote{Department of Economics,
University of Copenhagen; e-mail:
\tt{[email removed]}.}}
\maketitle
\begin{abstract}
We develop a new method for selecting the penalty parameter for $\ell_{1}$-penalized M-estimators in
high dimensions, which we refer to as bootstrapping after cross-validation. We derive rates of
convergence for the corresponding $\ell_1$-penalized M-estimator and also for the
post-$\ell_1$-penalized M-estimator, which refits the non-zero entries of the former estimator
without penalty in the criterion function. We demonstrate via simulations that our methods are not
dominated by cross-validation in terms of estimation errors and can outperform cross-validation in
terms of inference. As an empirical illustration, we revisit \citet{fryer_jr_empirical_2019}, who
investigated racial differences in police use of force, and confirm his findings.
\end{abstract}
\medskip
\noindent
\textbf{Keywords:} Penalty parameter selection, penalized M-estimation,
high-dimensional models, sparsity, cross-validation, bootstrap, inference,
one-step debiasing.
\section{Introduction}\label{sec:Introduction}
High-dimensional models have attracted substantial attention both in the econometrics and in the
statistics/machine learning literature, see e.g.~\citet{belloni2018highdimensional} and
\citet{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
\begin{equation}
\boldsymbol{\theta}_{0}=\operatornamewithlimits{argmin}\limits_{\boldsymbol{\theta}\in\Theta}\mathrm{E}[m(\boldsymbol{X}^{\top}\boldsymbol{\theta},\boldsymbol{Y})],\label{eq:EstimandIntro}
\end{equation}
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{eq:EstimandIntro}) also covers many other
cross-sectional models and associated modern as well as classical estimation approaches including
logit and probit models, logistic calibration \citep{tan2020regularized}, covariate balancing
\citep{imai_covariate_2014}, and expectile regression \citep{newey_asymmetric_1987}. It also
subsumes approaches to estimation of panel-data models such as the fixed-effects/conditional logit
for binary outcomes \citep{rasch1960probabilistic}, trimmed least-absolute-deviations and trimmed
least-squares for censored outcomes \citep{honore_trimmed_1992}, and partial likelihood approaches
to heterogeneous panel models for duration \citep{chamberlain_heterogeneity_1985}. We detail some of
these examples in Section \ref{sec:Examples}.
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{sec:Deterministic-Bounds},
approximate sparsity captures the idea that, even though the number of \emph{candidate} regressors
$p$ can be very large, the number of \emph{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)
\begin{equation}
\widehat{\boldsymbol{\theta}}\left(\lambda\right)\in\widehat{\Theta}\left(\lambda\right):=\operatornamewithlimits{argmin}\limits_{\boldsymbol{\theta}\in\Theta}\bigg\{\frac{1}{n}\sum_{i=1}^{n}m(\boldsymbol{X}_{i}^{\top}\boldsymbol{\theta},\boldsymbol{Y}_{i})+\lambda\Vert\boldsymbol{\theta}\Vert_{1}\bigg\},\label{eq:ell1PenalizedMEstimationIntro}
\end{equation}
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{sec:ExistenceSparsityAndUniqueness}.} 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 \eqref{eq:ell1PenalizedMEstimationIntro}.
Implementing the estimator $\widehat{\boldsymbol{\theta}}(\lambda)$ requires choosing $\lambda$. To do so, we
first extend a probabilistic bound from \citet{belloni_l1-penalized_2011}, obtained for
$\ell_{1}$-penalized quantile regression, to our general $\ell_{1}$-penalized M-estimation setting
(\ref{eq:ell1PenalizedMEstimationIntro}). (See also \citet{negahban_unified_2012} for independently
developed and closely related results.) The bound, which we state in Section
\ref{sec:Deterministic-Bounds}, 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
\begin{equation}
\lambda \geqslant c_{0}\max_{1\leqslant j\leqslant p}\bigg|\frac{1}{n}\sum_{i=1}^{n}m_{1}'(\boldsymbol{X}_{i}^{\top}\boldsymbol{\theta}_{0},\boldsymbol{Y}_{i})X_{i,j}\bigg|\label{eq:ScoreDominationIntro}
\end{equation}
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
\begin{equation}
q_n(1-\alpha):=\left(1-\alpha\right)\text{-quantile of }\max_{1\leqslant j\leqslant p}\bigg|\frac{1}{n}\sum_{i=1}^{n}m_{1}'(\boldsymbol{X}_{i}^{\top}\boldsymbol{\theta}_{0},\boldsymbol{Y}_{i})X_{i,j}\bigg|,\label{eq:ScoreQuantileIntro}
\end{equation}
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{eq:ScoreQuantileIntro}) 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 \citet{belloni2018highdimensional}, and to set $\lambda=c_{0}\widehat{q}(1-\alpha)$,
which we refer to as the \textit{bootstrap-after-cross-validation} (BCV) \textit{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{sec:Bootstrapping-the-Penalty}. In addition, we show in Section \ref{sec: inference} 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.~\citet{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 \citet{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{sec:
inference}.\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 \cite{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
\citet{chetverikov_cross-validated_2016} and \citet{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{sec:Simulations} for details.
Another alternative to our method is to base the penalty parameter choice on self-normalized
moderate deviation (SNMD) theory, as proposed in \cite{belloni_sparse_2012} for the linear mean
regression model and extended in \cite{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{sec:Examples}.
To showcase our method using real data, in Section \ref{sec:Application} we revisit the setting of
\citet{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 \citet{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 \citet{wainwright_high-dimensional_2019} and focus
here on only a few key references.
\citet{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{eq:ell1PenalizedMEstimationIntro}) 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, \citet[p.~621]{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.
\citet{negahban_unified_2012} develop error guarantees in a very general setting, and when
specialized to our setting (\ref{eq:ell1PenalizedMEstimationIntro}), their results become quite
similar to those in our Theorem \ref{thm:NonAsymptoticProbabilisticBounds}. The same authors also
note that a challenge to using these results in practice is that the random variable in
(\ref{eq:ScoreDominationIntro}) is usually impossible to compute because it depends on the unknown
vector $\boldsymbol{\theta}_{0}$ (\emph{ibid.,} p.~547). It is exactly this challenge that we overcome in this
paper. \citet{belloni_l1-penalized_2011} study the high-dimensional quantile regression model and
note that the distribution of the random variable in (\ref{eq:ScoreDominationIntro}) is in this case
pivotal, making the choice of the penalty parameter simple. Similarly, \cite{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{eq:ScoreDominationIntro}) is pivotal.\footnote{With a
known censoring propensity, the \cite{buchinsky1998alternative} linear programming estimator for
censored quantile regression boils down to a variant of quantile regression, thus leading to
pivotality.} Finally, \citet{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{sec:Examples} we
provide a portfolio of examples that constitute possible applications of our
method. In Section \ref{sec:Deterministic-Bounds} we develop bounds on the
estimation error of the $\ell_1$-ME, which motivate our
method for choosing the penalty parameter. In Section
\ref{sec:Bootstrapping-the-Penalty}, we introduce the BCV penalty method and
derive convergence rates for the resulting $\ell_1$-ME and
post-$\ell_1$-ME. In Section \ref{sec: inference}, we show how to perform
inference on individual components of $\boldsymbol{\theta}_0$ via debiasing. In Section
\ref{sec:Simulations}, 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{sec:Application}, we apply our method to the empirical
setting of \citet{fryer_jr_empirical_2019}. All proofs are relegated to the
Online Appendices.
\subsection*{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.
\section{Examples\label{sec:Examples}}
In this section, we discuss a variety of models that fit into the M-estimation
framework (\ref{eq:EstimandIntro}) 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{sec:Additional-Examples}.
\begin{example}
[\textbf{Binary Response Model}]\label{exa:Logit}
A relatively simple model fitting our framework is the \emph{binary} \emph{response model}, i.e. a
model for an outcome $Y\in\{0,1\} $ with
\[
\mathrm{P}(Y=1\mid\boldsymbol{X})=F(\boldsymbol{X}^{\top}\boldsymbol{\theta}_{0}),
\]
for a known cumulative distribution function (CDF) $F:\mathbb{R}\to\left(0,1\right)$. The log-likelihood of
this model yields the following loss function:
\begin{equation}
m\left(t,y\right)=-y\ln F\left(t\right)-\left(1-y\right)\ln\left(1-F\left(t\right)\right).\label{eq:LossBinaryResponse}
\end{equation}
The \emph{logit} model arises from setting
$F\left(t\right)=1/\left(1+\mathrm{e}^{-t}\right)=:\Lambda\left(t\right)$, the standard logistic
CDF, and the loss function reduces in this case to
\begin{equation}
m\left(t,y\right)=\ln\left(1+\mathrm{e}^{t}\right)-yt.\label{eq:LossLogit}
\end{equation}
The \emph{probit} model arises from setting
$F\left(t\right)=\int_{-\infty}^{t}\left(2\pi\right)^{-1/2}\allowbreak\mathrm{e}^{-u^{2}/2}\mathrm{d}u=:\Phi\left(t\right)$,
the standard normal CDF, and the loss function in this case becomes
\begin{equation}
m\left(t,y\right)=-y\ln\Phi\left(t\right)-\left(1-y\right)\ln\left(1-\Phi\left(t\right)\right).\label{eq:LossProbit}
\end{equation}
Both loss functions (\ref{eq:LossLogit}) and (\ref{eq:LossProbit}) are convex in $t$.
More generally, any binary response model with both $F$ and complementary CDF $1-F$ being
log-concave leads to a loss (\ref{eq:LossBinaryResponse}) that is convex in $t$. For these
log-concavities it suffices that $F$ admits a probability density function (PDF) $f=F'$, which is
itself positive and log-concave \citep[Section 5]{pratt_concavity_1981}. Both the standard logistic
and standard normal PDFs are log-concave. Also, $\ln f$ is concave whenever $f$ is of the (Subbotin)
form $f(t)\propto\mathrm{e}^{-\vert t\vert^{a}/a}$ for some $a\in[1,\infty)$, the extreme case being
the Laplace distribution. See \citet[Section 6]{pratt_concavity_1981} for additional examples. We
focus on the logit and probit cases for concreteness.\qed
\end{example}
\begin{example}[\textbf{Ordered Response Model}]\label{exa:LogConcaveOrderedResponse}
Consider the \emph{ordered response model}, i.e.~a model for an outcome $Y\in\{0,1,\dots,V\}$ with
\[
\mathrm{P}(Y=v\mid\boldsymbol{X})=F(\alpha_{v+1}-\boldsymbol{X}^{\top}\boldsymbol{\theta}_{0})-F(\alpha_{v}-\boldsymbol{X}^{\top}\boldsymbol{\theta}_{0}),\quad v\in\left[V\right] ,
\]
for a known CDF $F:\mathbb{R}\to(0,1)$ and known cut-off points
$-\infty=\alpha_{0}<\alpha_{1}<\cdots<\alpha_{V}<\alpha_{V+1}=+\infty$. (We interpret
$F\left(-\infty\right)$ as zero and $F\left(+\infty\right)$ as one to subsume the end cases.) The
log-likelihood of this model yields the loss function
\begin{equation}
m\left(t,y\right)=-\sum_{v=0}^{V}\mathbf{1}\left(y=v\right)\ln\left(F\left(\alpha_{v+1}-t\right)-F\left(\alpha_{v}-t\right)\right),\label{eq:LossOrderedResponse}
\end{equation}
which is convex in $t$ for any distribution $F$ admitting a positive and log-concave PDF $f=F'$
\citep[Section 3]{pratt_concavity_1981}. See Example \ref{exa:Logit} for specific distributions
satisfying this criterion. As for binary response, we focus on the logit and probit cases.\qed
\end{example}
\begin{example}
[\textbf{Expectile Model}]\label{exa:Expectile} \citet{newey_asymmetric_1987}
study the conditional ($\tau$th) \textit{expectile model}
$\mu_{\tau}(Y\mid\boldsymbol{X})=\boldsymbol{X}^{\top}\boldsymbol{\theta}_{0}$, where
$\tau\in(0,1)$ is a known number, and propose the\emph{ asymmetric least
squares} (ALS) estimator of $\boldsymbol{\theta}_{0}$ in this model. This estimator
can be understood as an M-estimator with loss of the form
\begin{equation}
m\left(t,y\right)=\rho_{\tau}\left(y-t\right),\label{eq: ALS}
\end{equation}
with $\rho_{\tau}:\mathbb{R}\to\mathbb{R}$ being the ``swoosh'' function given by
\[
\rho_{\tau}\left(u\right)=\left|\tau-\mathbf{1}\left(u<0\right)\right|u^{2}=\begin{cases}
\left(1-\tau\right)u^{2}, &\text{if}\; u<0,\\
\tau u^{2}, &\text{if}\; u\geqslant0,
\end{cases}
\]
a piecewise quadratic and continuously differentiable analogue of the ``check'' function known from
the quantile regression literature. The ALS estimator can be interpreted as a maximum likelihood
estimator when model disturbances arise from a normal distribution with unequal weights placed on
positive and negative disturbances \citep{aigner_estimation_1976, philipps2022mle}. Note that
$m(\cdot,y)$ in (\ref{eq: ALS}) is convex but not twice differentiable (at $y$) unless
$\tau=1/2$.\qed
\end{example}
\begin{example}
[\textbf{Panel Censored Model}]
\label{exa:PanelCensoredRegressionAndTrimming}
Consider the \textit{panel censored model}
\[
Y_{\tau}=\max\left(0,\gamma+\boldsymbol{X}_{\tau}^{\top}\boldsymbol{\theta}_0+\varepsilon_{\tau}\right),\quad \tau\in\{1,2\},
\]
where $\boldsymbol{Y}=(Y_{1},Y_{2})^{\top}\in[0,\infty)^2$ is a pair of outcome variables,
$(\boldsymbol{X}_{1}^{\top},\boldsymbol{X}_{2}^{\top})^{\top}$ is a vector of regressors, $\gamma$ is a unit-specific
(possibly random) unobserved fixed effect, and $\varepsilon_{1}$ and $\varepsilon_{2}$ are
unobserved error terms, which may or may not be centered. \citet{honore_trimmed_1992} shows that
under certain conditions, including exchangeability of $\varepsilon_{1}$ and $\varepsilon_{2}$
conditional on $(\boldsymbol{X}_{1},\boldsymbol{X}_{2},\gamma)$, $\boldsymbol{\theta}_{0}$ in this model can be identified by
$\boldsymbol{\theta}_{0}=\operatornamewithlimits{argmin}\limits\nolimits_{\boldsymbol{\theta}\in\mathbb{R}^{p}}\mathrm{E}[m(\boldsymbol{X}^{\top}\boldsymbol{\theta},\boldsymbol{Y})],$ with $\boldsymbol{X} := \boldsymbol{X}_1 -
\boldsymbol{X}_2$ and $m$ being the \emph{trimmed loss }function
\begin{equation}
m\left(t,\boldsymbol{y}\right)=\begin{cases}
\Xi\left(y_{1}\right)-\left(y_{2}+t\right)\xi\left(y_{1}\right), &\text{if}\; t\in\left(-\infty,-y_2\right],\\
\Xi\left(y_{1}-y_{2}-t\right), &\text{if}\; t\in\left(-y_{2},y_{1}\right),\\
\Xi\left(-y_{2}\right)-\left(t-y_{1}\right)\xi\left(-y_{2}\right), &\text{if}\; t\in\left[y_{1},\infty\right),
\end{cases}\label{eq:TrimmedLoss}
\end{equation}
and either $\Xi=\left\vert\cdot\right\vert$ or $\Xi=(\cdot)^{2}$ and $\xi$ its derivative (when
defined).\footnote{When $\Xi=\left|\cdot\right|$, we set $\xi\left(0\right):=0$ to make
\eqref{eq:TrimmedLoss} consistent with formulas in \citet{honore_trimmed_1992}.} These choices lead
to \emph{trimmed least absolute deviations} (trimmed LAD) and \emph{trimmed least squares} (trimmed
LS) estimators, respectively, both of which are based on loss functions convex in $t$. Note that
trimmed LAD is based on a non-differentiable loss $m(\cdot,\boldsymbol{y})$, while trimmed LS is based on a
continuously differentiable but not twice differentiable loss.\qed
\end{example}
\section{Non-Asymptotic Bounds on Estimation Error\label{sec:Deterministic-Bounds}}
In this section, we derive probabilistic bounds on the error of the
$\ell_{1}$-ME (\ref{eq:ell1PenalizedMEstimationIntro}) 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{as: diff and int}, \ref{assu:Margin}, and
\ref{assu:LossLocallyLipschitzAndMore} 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{sec:VerificationLinearModelSquareLoss}
and for all examples in Section \ref{sec:Examples} in Appendix \ref{sec:
verification}.
\begin{assumption}
[\textbf{Parameter Space}]\label{assu:ParameterSpace} The parameter space
$\Theta$ is a non-empty convex subset of $\mathbb{R}^{p}$ for which $\boldsymbol{\theta}_{0}$ is
interior.
\end{assumption}
\begin{assumption}
[\textbf{Convexity}]\label{assu:Convexity} The function
$m\left(\cdot,\boldsymbol{y}\right)$ is convex for all $\boldsymbol{y}\in\mathcal{Y}$.
\end{assumption}
\begin{assumption}
[\textbf{Differentiability and Integrability}]\label{as: diff and int} The
derivative $m_{1}'(\boldsymbol{X}^{\top}\boldsymbol{\theta},\boldsymbol{Y})$ exists almost surely for all
$\boldsymbol{\theta}\in\Theta$, and $\mathrm{E}[\lvertm(\boldsymbol{X}^{\top}\boldsymbol{\theta},\boldsymbol{Y})\rvert]<\infty$ for
all $\boldsymbol{\theta}\in\Theta$.
\end{assumption}
Assumption \ref{assu:ParameterSpace} 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{assu:Convexity} is satisfied in all examples from
the previous section, as discussed there. In the same examples, Assumption
\ref{as: diff and int} imposes minor integrability conditions on the random
vectors $\boldsymbol{X}$ and $\boldsymbol{Y}$. In addition, in the case of Example
\ref{exa:PanelCensoredRegressionAndTrimming} 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{sec: verification} for
details.
Further, define the \emph{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 \eqref{eq:EstimandIntro}, 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$.
\begin{assumption}
[\textbf{Margin}]\label{assu:Margin} There are constants $c_{M}\in(0,1]$ and $c_{M}'\in(0,\infty]$
such that for all $\boldsymbol{\theta}\in\Theta$ satisfying $\Vert\boldsymbol{\theta}-\boldsymbol{\theta}_{0}\Vert_{2}\leqslant c_{M}'$,
we have $\mathcal{E}\left(\boldsymbol{\theta}\right)\geqslant c_{M}\Vert\boldsymbol{\theta}-\boldsymbol{\theta}_{0}\Vert_{2}^{2}$.
\end{assumption}
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{exa:PanelCensoredRegressionAndTrimming}, it also requires $Y_1$ and $Y_2$ to be different with
positive probability. Also, our formal analysis reveals that Assumption \ref{assu:Margin} 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 \emph{sparse} vectors
$\boldsymbol{\theta}$. We have opted for a less general statement to avoid additional technicalities.
\begin{comment}
\begin{itemize}
\item J: New ($\ell_2$) margin nbh. Old ($\ell_1$) margin nbh was stronger than necessary.
\item J: Have allowed $c_{M}'=\infty$ to weaken side conditions in case
of square loss.
\end{itemize}
\end{comment}
The following assumption requires additional technical regularity of the loss function.
\begin{assumption}[\textbf{Local Loss}]\label{assu:LossLocallyLipschitzAndMore}
There are constants $c_{L}\in(0,\infty]$, $C_{L}\in[1,\infty)$ and $r\in(4,\infty)$, a non-random
sequence $B_n$ in $[1,\infty)$, and a function $L:\mathcal{X}\times\mathcal Y\to[1,\infty)$ such
that
\begin{enumerate}
\item \label{enu:LossLocallyLipschitz} for all
$(\boldsymbol{x},\boldsymbol{y})\in\mathcal{X}\times\mathcal Y$ and all $(t_1,t_2)\in\mathbb{R}^2$ satisfying
$\left\lvertt_1\right\rvert\lor\left\lvertt_2\right\rvert\leqslant c_{L},$
\begin{align}
\left\lvertm\left(\boldsymbol{x}^{\top}\boldsymbol{\theta}_{0}+t_1,\boldsymbol{y}\right)-m\left(\boldsymbol{x}^{\top}\boldsymbol{\theta}_{0}+t_2,\boldsymbol{y}\right)\right\rvert & \leqslant L\left(\boldsymbol{x},\boldsymbol{y}\right)\left\lvertt_1-t_2\right\rvert\label{eq:LocallyLipschitz}
\end{align}
with $\max_{1\leqslant j\leqslant p}\mathrm{E}[|L(\boldsymbol{X},\boldsymbol{Y})X_j|^2]\leqslant C_L^2$ and
$\mathrm{E}[|L(\boldsymbol{X},\boldsymbol{Y})\|\boldsymbol{X}\|_{\infty}|^r]\leqslant B_n^r$;
\item \label{enu:LossMeanSquareEll2Conts} for all $\boldsymbol{\theta}\in\Theta$ satisfying
$\|\boldsymbol{\theta} - \boldsymbol{\theta}_0\|_2\leqslant c_L$, we have
\begin{align*}
\mathrm{E}\left[\left|m\left(\boldsymbol{X}^{\top}\boldsymbol{\theta},\boldsymbol{Y}\right)-m\left(\boldsymbol{X}^{\top}\boldsymbol{\theta}_{0},\boldsymbol{Y}\right)\right|^2\right] & \leqslant C_L^2 \|\boldsymbol{\theta} - \boldsymbol{\theta}_0\|_2^2;
\end{align*}
\item \label{enu:ResidualMeanSquareEll2Conts} for all $\boldsymbol{\theta}\in\Theta$ satisfying $\|\boldsymbol{\theta} -
\boldsymbol{\theta}_0\|_2\leqslant c_L$, we have
\begin{equation}\label{eq: derivative mean square differentiability}
\mathrm{E}\left[\left| m_{1}'\left(\boldsymbol{X}^{\top}\boldsymbol{\theta},\boldsymbol{Y}\right)-m_{1}'\left(\boldsymbol{X}^{\top}\boldsymbol{\theta}_{0},\boldsymbol{Y}\right)\right|^{2}\right]\leqslant C_{L}^{2}\|\boldsymbol{\theta} - \boldsymbol{\theta}_0\|_2.
\end{equation}
\end{enumerate}
\end{assumption}
\begin{comment}
\begin{itemize}
\item J: $c_{L}=\infty$ allows global Lipschitzness (makes some sample size
qualifications vacuous).
\end{itemize}
\end{comment}
Assumption
\ref{assu:LossLocallyLipschitzAndMore}.\ref{enu:LossLocallyLipschitz} 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 \eqref{eq:LocallyLipschitz} actually
follows from the loss convexity in Assumption \ref{assu:Convexity} \citep[][Theorem
10.4]{rockafellar_convex_1970}, so Assumption
\ref{assu:LossLocallyLipschitzAndMore}.\ref{enu:LossLocallyLipschitz} should be regarded as a mild
moment condition. Assumptions
\ref{assu:LossLocallyLipschitzAndMore}.\ref{enu:LossMeanSquareEll2Conts} and
\ref{assu:LossLocallyLipschitzAndMore}.\ref{enu:ResidualMeanSquareEll2Conts} 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 \emph{globally} Lipschitz uniformly in $\boldsymbol{y}$ (thus
allowing the choice $c_L=\infty$), Assumption
\ref{assu:LossLocallyLipschitzAndMore}.\ref{enu:LossLocallyLipschitz} boils down to the regressors
having sufficiently many absolute moments, and Assumption
\ref{assu:LossLocallyLipschitzAndMore}.\ref{enu:LossMeanSquareEll2Conts} 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.~\citet[Condition A.2]{belloni_new_2015} and \citet[Assumption 5]{soerensen_2024}.} Examples of
globally Lipschitz losses are the logit likelihood loss in Example \ref{exa:Logit} and the trimmed
LAD loss in Example \ref{exa:PanelCensoredRegressionAndTrimming}.
\begin{assumption}
[\textbf{Approximate Sparsity}]\label{assu:Approximate-Sparsity}
There is a constant $q\in\left[0,1\right]$ and a non-random sequence $s_{q}:=s_{q,n}$ in
$[1,\infty)$ such that $\sum_{j=1}^{p}|\theta_{0,j}|^{q}\leqslant s_{q}.$
\end{assumption}
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 \textit{exact} sparsity with (at most) $s_{0}$
non-zero entries. When $q>0$, we have only \textit{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{rem:SparsityNotions} for further discussion.
Under Assumption \ref{as: diff and int}, we can (almost surely) define
$\boldsymbol{S}_n\in\mathbb{R}^p$ by
\begin{equation}
\boldsymbol{S}_n:=\mathbb{E}_{n}\left[\left.\frac{\partial}{\partial\boldsymbol{\theta}}m(\boldsymbol{X}_{i}^{\top}\boldsymbol{\theta},\boldsymbol{Y}_{i})\right|_{\boldsymbol{\theta} = \boldsymbol{\theta}_0}\right]=\mathbb{E}_{n}\left[m_1'(\boldsymbol{X}_{i}^{\top}\boldsymbol{\theta}_0,\boldsymbol{Y}_{i})\boldsymbol{X}_{i}\right].\label{eq:Score}
\end{equation}
In this paper we refer to $\boldsymbol{S}_n$ as the \emph{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{sec:ProofsDeterministicBounds}, builds on arguments
of \cite{belloni_l1-penalized_2011}. Related statements appear also in
\cite{van_de_geer_high-dimensional_2008}, \cite{bickel_simultaneous_2009}, and
\cite{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}$.
\begin{thm}
[\textbf{Non-Asymptotic Error Bounds for $\boldsymbol{\ell_1}$-ME}]\label{thm:NonAsymptoticProbabilisticBounds} Let Assumptions
\ref{assu:ParameterSpace}--\ref{assu:Approximate-Sparsity} hold, let
$\overline{\lambda}_{n}$ be a non-random sequence in $(0,\infty)$, let $c_0\in(1,\infty)$, and define
\begin{align*}
\ u_{n} & :=\frac{4c_{0}\sqrt{s_{q}\eta_{n}^{-q}}}{(c_{0}-1)c_M}\left(C_{L}\eta_{n}+\overline{\lambda}_{n}\right).
\end{align*}
Then there is a universal constant $C\in[1,\infty)$ such that for all $n\in \mathbb{N}$
and $t\in[1,\infty)$ satisfying
$$
\eta_{n}\leqslant1,\
Cu_{n}\leqslant c_{M}', \
\frac{B_n^2\ln(pn)}{\sqrt n}\leqslant C_L^2 \text{ and } \
t n^{1/r} B_n\left(Cu_{n}\sqrt{s_{q}\eta_{n}^{-q}}+s_{q}\eta_{n}^{1-q}\right)\leqslant\frac{\left(c_{0}-1\right)c_{L}}{2c_{0}},
$$
we have
\begin{align*}
\sup_{\mathclap{\widehat{\boldsymbol{\theta}}\in\widehat{\Theta}\left(\lambda\right)}}\|\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}\|_{2} & \leqslant Cu_{n}
\quad\text{and}\quad\sup_{\mathclap{\widehat{\boldsymbol{\theta}}\in\widehat{\Theta}\left(\lambda\right)}}\|\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}\|_{1}\leqslant\frac{2c_{0}}{c_{0}-1}\left(Cu_{n}\sqrt{s_{q}\eta_{n}^{-q}}+s_{q}\eta_{n}^{1-q}\right)
\end{align*}
with probability at least
$1-\mathrm{P}(\lambda<c_{0}\|\boldsymbol{S}_n\|_{\infty})-\mathrm{P}(\lambda>\overline{\lambda}_{n})-
4t^{-r} - C/\ln^2(pn) - n^{-1}.$
\end{thm}
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{thm:NonAsymptoticProbabilisticBounds} for the post-$\ell_1$-ME
in Appendix \ref{sec:analysis post estimator}, but the general principle for
choosing $\lambda$ remains the same.
\begin{rem}[\textbf{Non-Uniqueness}]
Like similar statements appearing in the literature, Theorem
\ref{thm:NonAsymptoticProbabilisticBounds} concerns the entire set $\widehat\Theta(\lambda)$ of
optimizers for the convex minimization problem (\ref{eq:ell1PenalizedMEstimationIntro}). While the
objective function is presumed convex, it need not be strictly convex, and the global minimum may be
attained at more than one point. For example, no matter the choice of $\Xi$, the trimmed loss
function \eqref{eq:TrimmedLoss} in Example \ref{exa:PanelCensoredRegressionAndTrimming} will have
linear pieces and need therefore not produce a strictly convex objective function. The bounds stated
here (and in what follows) hold for any of these optimizers. See also Appendix
\ref{sec:ExistenceSparsityAndUniqueness} for sufficient conditions for solution existence and
uniqueness as well as related (sparsity) properties. Despite the possible multiplicity, we sometimes
refer to any element $\widehat{\boldsymbol{\theta}}\in\widehat{\Theta}\left(\lambda\right)$ as \emph{the}
$\ell_1$-ME. \qed
\end{rem}
\begin{rem}[\textbf{Margin}]
Our convexity, interiority and differentiability assumptions suffice to show
that the excess risk function $\mathcal E(\boldsymbol{\theta})$ is differentiable at
$\boldsymbol{\theta}_0$, and so the estimand $\boldsymbol{\theta}_{0}$ must satisfy the population
first-order condition $\nabla \mathcal E\left(\boldsymbol{\theta}_{0}\right)=\mathbf{0}$.
Assumption \ref{assu:Margin} therefore amounts to assuming that $\mathcal
E(\boldsymbol{\theta})$ admits a quadratic \emph{margin} near $\boldsymbol{\theta}_{0}$. The name
\emph{margin condition} appears to originate from \citet[Assumption
A1]{tsybakov_optimal_2004}, who invokes a similar assumption in a classification
context. \citet[Assumption B]{van_de_geer_high-dimensional_2008} contains a more
general formulation of margin behavior for estimation purposes. We consider the
(focal) quadratic case for the sake of simplicity.\qed
\end{rem}
\begin{rem}
[\textbf{Sparsity Notions}]\label{rem:SparsityNotions} In \citet[Section
4.3]{negahban_unified_2012} the sparsity in Assumption
\ref{assu:Approximate-Sparsity} is referred to as \emph{strong} for $q=0$ and
\emph{weak }for $q>0.$ \citet[Chapter 7]{wainwright_high-dimensional_2019}
distinguishes between \emph{strong }$\ell_{q}$-balls (like
$\{\boldsymbol{\theta}\in\mathbb{R}^{p};\sum_{j=1}^{p}|\theta_{j}|^{q}\leqslant s_{q}\}$
implicitly considered here) and \emph{weak }$\ell_{q}$-balls, which impose a
polynomial decay in the non-increasing rearrangement of the absolute values of
the coefficients. In \citet{belloni2018highdimensional}, restricting $\boldsymbol{\theta}_{0}$ to a weak $\ell_{q}$-ball is referred to as \emph{approximate
sparsity}, and a $\boldsymbol{\theta}_{0}$ having bounded $\ell_{1}$ norm (i.e.~belonging to
a strong $\ell_{1}$-ball) is called \emph{dense}. Both strong $(q>0)$ and weak
ball restrictions formalize the idea of ``weak'' or ``approximate''
sparsity.\qed
\end{rem}
\begin{rem}[\textbf{Free Parameter}]
The free parameter $c_0\in(1,\infty)$ in Theorem \ref{thm:NonAsymptoticProbabilisticBounds} serves
as a trade-off between the likelihood of score domination on the one hand and the bound quality on
the other. A smaller $c_0\in(1,\infty)$ makes the event $\lambda\geqslant c_{0}\|\boldsymbol{S}_n\|_{\infty}$
more probable but also worsens the bounds. Note that the free parameter $c_0$ appears, either
explicitly or implicitly, in existing bounds as well.\footnote{A free parameter is explicit in both
\cite{belloni_l1-penalized_2011} and \cite{van_de_geer_high-dimensional_2008}. In deriving their
bounds both \cite{bickel_simultaneous_2009} (for the LASSO) and \cite{negahban_unified_2012} set
$c_0=2$.} While asymptotic theory provides no guidance on the choice of $c_0$, our
finite-sample experiments in Section \ref{sec:Simulations} indicate that increasing $c_0$ away from
one worsens performance but setting $c_0$ to any value near one, including one itself, does not
impact the results by much (cf.~Figures \ref{fig:MeanEll2ErrorBCVVaryingc0BCCHrule} and
\ref{fig:MeanEll2ErrorPostBCVVaryingc0BCCHrule}). Similar observations were made by \citet[Footnote
7]{belloni_sparse_2012} in the context of the LASSO.\qed
\end{rem}
\section{Bootstrapping after Cross-Validation\label{sec:Bootstrapping-the-Penalty}}
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{sub: boot penalty level}, 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{sec:ResidualEstimationViaCrossValidation}, 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
\emph{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.}
\subsection{Bootstrapping the Penalty Level}\label{sub: boot 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{as: diff and
int}. 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 \citet{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
\begin{equation}
\sup_{A\in\mathcal{A}_{p}}\left|\mathrm{P}\left(\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\boldsymbol{Z}_i\in A\right)-\mathrm{P}(\boldsymbol{\mathcal{N}}_n\in A)\right|\leqslant C_{b}\left(\frac{\widetilde{B}_n^4\ln^{7}\left(pn\right)}{n}\right)^{1/6},\label{eq:CCK_HDCLT}
\end{equation}
and, with probability approaching one,
\begin{equation}
\sup_{A\in\mathcal{A}_{p}}\left|\mathrm{P}\left(\left.\frac{1}{\sqrt{n}}\sum_{i=1}^{n}e_{i}\boldsymbol{Z}_{i}\in A\right|\left\{\boldsymbol{Z}_{i}\right\} _{i=1}^{n}\right)-\mathrm{P}(\boldsymbol{\mathcal{N}}_n\in A)\right|\leqslant C_{b}\left(\frac{\widetilde{B}_n^4\ln^{7}\left(pn\right)}{n}\right)^{1/6},\label{eq:CCK_HDBootstrap}
\end{equation}
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{eq:CCK_HDBootstrap}) 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 \eqref{eq: residual estimation side condition}. Suppose
therefore that residual estimators $\{\widehat{U}_{i}\}_{i=1}^{n}$ are
available. We then compute
\begin{equation}
\widehat{q}^{\texttt{bm}}\left(1-\alpha\right):=\left(1-\alpha\right)\text{-quantile of }\max_{1\leqslant j\leqslant p}\big|\mathbb{E}_{n}\big[e_{i}\widehat{U}_{i}X_{i,j}\big]\big|\text{ given }\{(\boldsymbol{X}_{i},\boldsymbol{Y}_i,\widehat{U}_{i})\}_{i=1}^{n},\label{eq:qhat1minusAlpha}
\end{equation}
and a penalty level follows as
\begin{equation}
\widehat{\lambda}^{\mathtt{bm}}_{\alpha}:=c_{0}\widehat{q}^{\texttt{bm}}\left(1-\alpha\right).\label{eq:BootstrapPenaltyLevel}
\end{equation}
We refer to this method for obtaining a penalty level as the \emph{bootstrap
method} (BM) and to $\widehat{\lambda}^{\mathtt{bm}}_{\alpha}$ itself as the
\emph{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.
\begin{assumption}
[\textbf{Residuals}]\label{assu:ResidualBootstrapMethod}
There are constants $c_{U}\in(0,\infty)$ and $C_U\in[1,\infty)$ and a non-random
sequence $\widetilde{B}_n$ in $[1,\infty)$, such that
\begin{inparaenum}[(1)]
\item\label{enu:ZijSecondMomentsBndAwayZero}
$c_U^2\leqslant\mathrm{E}[|UX_j|^2]\leqslant C_{U}^2$ for all $j\in[p]$,
\item\label{enu:Zj2pluskMomentBnd} $\mathrm{E}[|UX_j|^{4}]\leqslant\widetilde{B}_n^2$
for all $j\in[p]$, and
\item\label{enu:maxZijFourthMomentBnd}
$\mathrm{E}\left[\|U\boldsymbol{X}\|_{\infty}^4\right]\leqslant \widetilde{B}_n^4.$\end{inparaenum}
\end{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{sec: verification} for the examples in Section
\ref{sec:Examples}.
Our next result provides convergence rates for the $\ell_1$-ME based on the bootstrap method.
\begin{lem}
[\textbf{Convergence Rates: Generic Bootstrap Method}]\label{thm:RatesBM} Let
Assumptions \ref{assu:ParameterSpace}--\ref{assu:Approximate-Sparsity} and
\ref{assu:ResidualBootstrapMethod} hold, let $\delta_n$ be a non-random sequence
in $[0,\infty)$ such that
\begin{equation}\label{eq: residual estimation side condition}
\mathrm{P}\Big(\mathbb{E}_{n}[(\widehat{U}_{i}-U_{i})^{2}] > \delta_{n}^{2}/\ln^{2}\left(pn\right)\Big)\to 0,
\end{equation}
let $\widehat{\Theta}(\widehat{\lambda}^{\mathtt{bm}}_{\alpha})$ be the
solutions to the $\ell_{1}$-penalized M-estimation problem
(\ref{eq:ell1PenalizedMEstimationIntro}) with penalty level
$\lambda=\widehat{\lambda}^{\mathtt{bm}}_{\alpha}$ given in
(\ref{eq:BootstrapPenaltyLevel}) and $\alpha=\alpha_n\in(0,1)$ satisfying $\alpha_n\to0$
and $\ln(1/\alpha_n)\lesssim \ln(pn)$, and suppose that
\begin{equation}\label{eq: simple restriction bm}
n^{1/r}B_{n}\left(\delta_{n} + s_{q}\eta_{n}^{1-q}\right)\to0,\quad
\frac{B_{n}^{4}\ln^2(pn)}{n}\to0\quad\text{and}\quad
\frac{\widetilde{B}_{n}^{4}\ln^{7}\left(pn\right)}{n}\to0.
\end{equation}
Then
\[
\sup_{\mathclap{\widehat\boldsymbol{\theta}\in\widehat{\Theta}(\widehat{\lambda}^{\mathtt{bm}}_{\alpha})}}\|\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}\|_{2} \lesssim_{\mathrm{P}} \sqrt{s_{q}\eta_{n}^{2-q}}\quad\text{and}\quad
\sup_{\mathclap{\widehat\boldsymbol{\theta}\in\widehat{\Theta}(\widehat{\lambda}^{\mathtt{bm}}_{\alpha})}}\|\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}\|_{1}\lesssim_{\mathrm{P}} s_{q}\eta_{n}^{1-q}.
\]
\end{lem}
The idea of using a bootstrap procedure to select the penalty level in
high-dimensional estimation is in itself not new.
\cite{chernozhukov_gaussian_2013} use a Gaussian multiplier bootstrap to tune
the Dantzig selector \citep{candes_dantzig_2007} for the high-dimensional linear
model allowing both non-Gaussian and heteroskedastic errors. Note, however, that
\citet[Theorem 4.2]{chernozhukov_gaussian_2013} presumes access to a preliminary
Dantzig selector, which is used to estimate residuals. The condition (\ref{eq:
residual estimation side condition}) 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{rem:ConvergenceRates} and
\ref{rem:DenseCase}.
\subsection{Cross-Validating Residuals\label{sec:ResidualEstimationViaCrossValidation}}
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
\begin{equation}
I_{k}=\left\{ \left(k-1\right)n/K+1,\dotsc,kn/K\right\} ,\quad k\in\left[K\right] ,\label{eq:CVEvenPartition}
\end{equation}
is natural, but not necessary. For the formal results below, we only require
that each $I_{k}$ specifies a ``substantial'' subsample; see Assumption
\ref{assu:DataPartition} below.
Let $\Lambda_{n}$ denote a finite subset of $(0,\infty)$ composed by candidate
penalty levels. In Assumption \ref{assu:CandidatePenalties} 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
\begin{equation}\label{eq:SubsampleEstimator}
\widehat{\boldsymbol{\theta}}_{I_{k}^{c}}\left(\lambda\right) \in \widehat{\Theta}_{I_{k}^{c}}\left(\lambda\right) := \operatornamewithlimits{argmin}\limits_{\boldsymbol{\theta}\in\Theta}\left\{ \mathbb{E}_{I_k^c}[m(\boldsymbol{X}_{i}^{\top}\boldsymbol{\theta},\boldsymbol{Y}_{i})]+\lambda\left\lVert\boldsymbol{\theta}\right\rVert_{1}\right\},
\end{equation}
for each candidate penalty level $\lambda\in\Lambda_{n}$ and holding out each
subsample $k\in[K]$ in turn. Second, determine the penalty level
\begin{equation}\label{eq:LambdaCVDefn}
\widehat{\lambda}^{\mathtt{cv}}\in\operatornamewithlimits{argmin}\limits_{\lambda\in\Lambda_{n}}\sum_{k=1}^{K}\sum_{i\in I_{k}}m\big(\boldsymbol{X}_{i}^{\top}\widehat{\boldsymbol{\theta}}_{I_{k}^{c}}\left(\lambda\right),\boldsymbol{Y}_{i}\big)
\end{equation}
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.,
\begin{equation}\label{eq:ResidualEstimatorCV}
\widehat{U}_{i}^{\mathtt{cv}}:=m_{1}'\big(\boldsymbol{X}_{i}^{\top}\widehat{\boldsymbol{\theta}}_{I_{k}^{c}}(\widehat{\lambda}^{\mathtt{cv}}),\boldsymbol{Y}_{i}\big),\quad i\in I_{k},\quad k\in\left[K\right].
\end{equation}
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{as: diff and int}. 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) \emph{method} for estimating
the quantile $q_n(1-\alpha)$,
\begin{equation}
\widehat{q}^{\mathtt{bcv}}\left(1-\alpha\right):=\left(1-\alpha\right)\text{-quantile of }\max_{1\leqslant j\leqslant p}\big|\mathbb{E}_{n}\big[e_{i}\widehat{U}_{i}^{\mathtt{cv}}X_{i,j}\big]\big|\text{ given }\{(\boldsymbol{X}_{i},\boldsymbol{Y}_i)\}_{i=1}^{n},\label{eq:qhat1minusAlpha CV}
\end{equation}
and the \emph{BCV penalty level} follows as
\begin{equation}
\widehat{\lambda}^{\mathtt{bcv}}_{\alpha}:=c_{0}\widehat{q}^{\mathtt{bcv}}\left(1-\alpha\right).\label{eq:BootstrapPenaltyLevel CV}
\end{equation}
To analyze the $\ell_1$-ME implied by BCV, we invoke the following two assumptions.
\begin{assumption}
[\textbf{Data Partition}]\label{assu:DataPartition} The number of folds $K\in\left\{
2,3,\dotsc\right\}$ is constant. There is a constant
$c_{D}\in\left(0,1\right)$ such that $\min_{1\leqslant k\leqslant
K}\left|I_{k}\right|\geqslant c_{D}n$.
\end{assumption}
\begin{assumption}
[\textbf{Candidate Penalties}]\label{assu:CandidatePenalties} There are
constants $c_{\Lambda}$ and $C_{\Lambda}$ in $(0,\infty)$ and
$a\in\left(0,1\right)$ such that
\[
\Lambda_{n}=\left\{ C_{\Lambda}a^{\ell};a^{\ell}\geqslant c_{\Lambda}/n,\ell\in\left\{ 0,1,2,\dotsc\right\} \right\} .
\]
\end{assumption}
Assumption \ref{assu:DataPartition} 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{assu:CandidatePenalties} 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{cor: convergence rate bootstrap after cv}
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
\eqref{eq: residual estimation side condition}. Combining this finding with Lemma \ref{thm:RatesBM},
we obtain convergence rates for the $\ell_1$-ME implied by BCV.
\begin{thm}
[\textbf{Convergence Rates: BCV Method, Penalized Estimator}]\label{cor: convergence rate bootstrap after cv} Let Assumptions
\ref{assu:ParameterSpace}--\ref{assu:Approximate-Sparsity} and
\ref{assu:ResidualBootstrapMethod}--\ref{assu:CandidatePenalties} hold, let
$\widehat\Theta(\widehat \lambda_{\alpha}^{\mathtt{bcv}})$ be the
solutions to the $\ell_1$-penalized M-estimation problem
\eqref{eq:ell1PenalizedMEstimationIntro} with penalty level $\lambda =
\widehat\lambda_{\alpha}^{\mathtt{bcv}}$ given in
\eqref{eq:BootstrapPenaltyLevel CV} and $\alpha=\alpha_n\in(0,1)$ satisfying $\alpha_n\to0$ and
$\ln(1/\alpha_n)\lesssim \ln(pn)$, and suppose that
\begin{equation}\label{eq: difficult restriction}
n^{1/r}B_ns_q\eta_n^{1-q}\to0,\quad
\frac{B_n^4 s_q(\ln(pn))^{5-q/2}(\ln n)^2}{n^{1-q/2-4/r}}\to0\quad\text{and}\quad
\frac{\widetilde{B}_n^4\ln^7\left(pn\right)}{n}\to0.
\end{equation}
Then
\begin{equation}\label{eq: rates ell1 me}
\sup_{\mathclap{\widehat\boldsymbol{\theta}\in\widehat{\Theta}(\widehat{\lambda}^{\mathtt{bcv}}_{\alpha})}}\|\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}\|_{2}\lesssim_{\mathrm{P}} \sqrt{s_{q}\eta_{n}^{2-q}}\quad\text{and}\quad
\sup_{\mathclap{\widehat\boldsymbol{\theta}\in\widehat{\Theta}(\widehat{\lambda}^{\mathtt{bcv}}_{\alpha})}}\|\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}\|_{1}\lesssim_{\mathrm{P}} s_{q}\eta_{n}^{1-q}.
\end{equation}
\end{thm}
\begin{rem}[\textbf{Convergence Rates}]\label{rem:ConvergenceRates}
The Theorem \ref{cor: convergence rate bootstrap after cv} (and Lemma \ref{thm:RatesBM}) convergence
rates are as one would expect in high-dimensional settings. For example, the $\ell_2$ rate in
\eqref{eq: rates ell1 me} coincides with that obtained for the LASSO in \citet[Corollary
3]{negahban_unified_2012} in the context of linear mean regression with $\boldsymbol{\theta}_0$ belonging to an
$\ell_q$-ball (Assumption \ref{assu:Approximate-Sparsity}). The rate is known to be minimax optimal
in the context of sparse linear mean regression \citep{ye2010rate,raskutti12}. We expect it to
remain optimal in the general high-dimensional M-estimation framework as well. In the special case
of exact sparsity $(q=0)$, the $\ell_2$ and $\ell_1$ rates in \eqref{eq: rates ell1 me} become
$\sqrt{s_0\ln(pn)/n}$ and $\sqrt{s_0^2\ln(pn)/n}$, respectively.\qed
\end{rem}
\begin{rem}[{\textbf{Dense Case}}]\label{rem:DenseCase}
In the dense case $(q=1)$ $\eta_n^{1-q}$ does not vanish, and so the side condition $n^{1/r}B_n s_q
\eta_n^{1-q}\to0$ in \eqref{eq: difficult restriction} fails. However, inspection of the proof
reveals that we actually require $n^{1/r}B_n s_q \eta_n^{1-q} / c_L\to0$ for $c_L\in(0,\infty]$ in
Assumption \ref{assu:LossLocallyLipschitzAndMore}. The latter condition is trivially satisfied when
$c_L=\infty$, which is allowed when the loss $m(t,\boldsymbol{y})$ is globally Lipschitz in its first argument
uniformly in $\boldsymbol{y}\in\mathcal Y$. Hence, provided the loss is globally Lipschitz, even in the dense
case Theorem \ref{cor: convergence rate bootstrap after cv} can produce the $\ell_2$ rate of
convergence $s_1^{1/2}(\ln(pn)/n)^{1/4}$. Examples of globally Lipschitz losses are the logit
likelihood loss in Example \ref{exa:Logit} and the trimmed LAD loss in Example
\ref{exa:PanelCensoredRegressionAndTrimming}. The side condition $n^{1/r}B_n s_q \eta_n^{1-q}\to0$
may also be relaxed in the special case of generalized linear models; see \citet[Section
4.4]{negahban_unified_2012} and \citet[Section 9.5]{wainwright_high-dimensional_2019} for
details.\qed
\end{rem}
\begin{rem}[\textbf{Model Sparsity and Regressor Regularity}]
The side condition
\[
\frac{B_n^4 s_q(\ln(pn))^{5-q/2}(\ln n)^2}{n^{1-q/2-4/r}}\to 0
\]
in \eqref{eq: difficult restriction} necessitates $r>8/(2-q)$, which reveals an
interplay between the model sparsity as captured by the constant $q$ in
Assumption \ref{assu:Approximate-Sparsity} and the regressor integrability as
captured by the constant $r$ in Assumption
\ref{assu:LossLocallyLipschitzAndMore}.\ref{enu:LossLocallyLipschitz}. In the
special case of exact sparsity $(q=0)$ the regressors are required to have more
than four finite moments.\qed
\end{rem}
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 \eqref{eq:ell1PenalizedMEstimationIntro}, 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
\begin{equation}\label{eq:post-ell1-ME-thetabar}
\widetilde\Theta(\mathrm{supp}(\overline\boldsymbol{\theta}))
:= \operatornamewithlimits{argmin}\limits_{\mathclap{\substack{\boldsymbol{\theta}\in\Theta,\\\supp\boldsymbol{\theta}\subseteq\mathrm{supp}(\overline\boldsymbol{\theta})}}}\mathbb{E}_{n}[m(\boldsymbol{X}_i^\top\boldsymbol{\theta},\boldsymbol{Y}_i)].
\end{equation}
Then for any $\ell_1$-ME, i.e.~a solution $\widehat\boldsymbol{\theta}\in\widehat\Theta(\lambda)$ to the
optimization problem in \eqref{eq:ell1PenalizedMEstimationIntro}, 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
\emph{all} post-$\ell_1$-MEs, which we denote
\begin{equation}\label{eq:AllPostEstimators}
\widetilde\Theta(\lambda):=
\bigcup_{\mathclap{\widehat\boldsymbol{\theta}\in\widehat\Theta(\lambda)}} \widetilde\Theta(\mathrm{supp}(\widehat\boldsymbol{\theta})).
\end{equation}
To analyze the post-$\ell_1$-ME, we will use the following two additional assumptions.
\begin{assumption}[\textbf{Smoothness}]\label{as: post estimator smoothness} The function
$m(\cdot,\boldsymbol{y})$ is differentiable for all $\boldsymbol{y}\in\mathcal Y$ with its derivative being
Lipschitz-continuous, i.e.~$|m_{1}'(t_2,\boldsymbol{y}) - m_1'(t_1,\boldsymbol{y})|\leqslant C_m|t_2 - t_1|$ for all
$(t_1,t_2,\boldsymbol{y})\in\mathbb R\times\mathbb R\times\mathcal Y$ and some constant $C_m\in(0,\infty)$.
\end{assumption}
\begin{assumption}[\textbf{Moments}]\label{as: post estimator moments}
There is a constant
$C_{ev}\in[1,\infty)$ such that $\mathrm{E}[(\boldsymbol{X}^\top \boldsymbol{\delta})^4]\leqslant
C_{ev}\|\boldsymbol{\delta}\|_2^4$ for all $\boldsymbol{\delta}\in\mathbb R^p$.
\end{assumption}
Assumption \ref{as: post estimator smoothness} strengthens the almost-sure differentiability in
Assumption \ref{as: diff and int}. The stronger smoothness requirement precludes the trimmed LAD
loss function in Example \ref{exa:PanelCensoredRegressionAndTrimming}, 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{sec:Examples}; see Appendix \ref{sec: verification}. Assumption \ref{as: post
estimator moments} 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.
\begin{thm}[\textbf{Convergence Rates: BCV Method, Post-Penalized Estimator}]
\label{thm: convergence rates post estimator}
Let Assumptions
\ref{assu:ParameterSpace}--\ref{assu:Approximate-Sparsity} and
\ref{assu:ResidualBootstrapMethod}--\ref{as: post estimator moments} hold, let
$\widetilde\Theta(\widehat \lambda_{\alpha}^{\mathtt{bcv}})$ be the set of post-$\ell_1$-penalized M-estimators \eqref{eq:AllPostEstimators} with penalty level $\lambda =
\widehat\lambda_{\alpha}^{\mathtt{bcv}}$ given in
\eqref{eq:BootstrapPenaltyLevel CV} and $\alpha=\alpha_n\in(0,1)$ satisfying $\alpha_n\to0$ and
$\ln(1/\alpha_n)\lesssim \ln(pn)$, and suppose that
\begin{equation}\label{eq: even more difficult restriction}
n^{1/r}B_ns_q\eta_n^{1-q}\ln(p n)\to0,\ \frac{B_n^4 s_q(\ln(pn))^{5-q/2}(\ln n)^2}{n^{1-q/2-4/r}}\to 0\ \text{and } \;\frac{\widetilde{B}_n^4\ln^7\left(pn\right)}{n}\to 0.
\end{equation}
Then
\begin{equation}\label{eq: rates post ell1 me}
\sup_{\mathclap{\widetilde\boldsymbol{\theta}\in\widetilde{\Theta}(\widehat{\lambda}^{\mathtt{bcv}}_{\alpha})}}\|\widetilde{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}\|_{2}\lesssim_{\mathrm{P}} \sqrt{s_{q}\eta_{n}^{2-q}\ln(pn)}\quad\text{and}\quad
\sup_{\mathclap{\widetilde\boldsymbol{\theta}\in\widetilde{\Theta}(\widehat{\lambda}^{\mathtt{bcv}}_{\alpha})}}\|\widetilde{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}\|_{1}\lesssim_{\mathrm{P}} s_{q}\eta_{n}^{1-q}\ln(pn).
\end{equation}
\end{thm}
\begin{rem}[{\textbf{Comparison of Rates for $\ell_1$-ME and Post-$\ell_1$-ME}}]
The convergence rates for the post-$\ell_1$-ME we derive here are slightly slower than those we
derived for the $\ell_1$-ME itself in Theorem \ref{cor: convergence rate bootstrap after cv}. The
technical reason for this difference is that the analysis of the post-$\ell_1$-ME requires not only
that the penalty parameter $\lambda$ is not too \textit{large} but also that it is not too
\textit{small}, as we may end up with ``too many'' selected variables; see Appendix
\ref{sec:analysis post estimator} for details. In turn, our BCV method may yield low values of the
penalty level if there is substantial correlation between regressors in the vector $\boldsymbol{X}$. A simple
solution to this issue would be to censor the BCV penalty parameter
$\widehat\lambda_{\alpha}^{\mathtt{bcv}}$ from below so that it shrinks to zero no faster than
$\eta_n$, i.e.~to replace $\widehat\lambda_{\alpha}^{\mathtt{bcv}}$ by
$\max\{\widehat\lambda_{\alpha}^{\mathtt{bcv}}, c\eta_n\}$ for some small user-chosen constant
$c\in(0,\infty)$. In this case, the rates in \eqref{eq: rates post ell1 me} would coincide with the
rates in \eqref{eq: rates ell1 me}. However, we prefer to state a slightly slower rate, as in
\eqref{eq: rates post ell1 me}, over the necessity to introduce another tuning parameter $c$ (for
which there is no obvious guiding principle). Moreover, under the additional assumption that the
elements of the regressor vector $\boldsymbol{X}$ are not too correlated, it is possible to derive the same
rates as in \eqref{eq: rates ell1 me} for the post-$\ell_1$-ME even without censoring, in which case
the rates for post-$\ell_1$-ME are as good as those for $\ell_1$-ME.\qed
\end{rem}
\section{Debiased Estimation and Inference}\label{sec: 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 \eqref{eq:EstimandIntro}.
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
\citep{chernozhukov2018double}. We note that the tools developed in this section rule out the
trimmed LAD loss in Example \ref{exa:PanelCensoredRegressionAndTrimming}, as this function is not
sufficiently smooth to satisfy our Assumption \ref{as: smoothness inference}.
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:
\begin{equation}\label{eq: definition of mu0}
\mathrm{E}[m''_{11}(\boldsymbol{X}^{\top}\boldsymbol{\theta}_0,\boldsymbol{Y})(D-\boldsymbol{W}^{\top}\boldsymbol{\mu}_0)\boldsymbol{W}] = \mathbf{0}_{p-1}.
\end{equation}
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{lem:ExistenceAndUniquenessOfMu0} in the appendix for the precise
statement.} With this definition in mind, by inspecting the first-order
conditions associated with \eqref{eq:EstimandIntro}, we have
\begin{equation}\label{eq: estimating equation}
\mathrm{E}[m'_1(D\beta_0 + \boldsymbol{W}^{\top}\boldsymbol{\gamma}_0,\boldsymbol{Y})(D-\boldsymbol{W}^{\top}\boldsymbol{\mu}_0)] = 0.
\end{equation}
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 \eqref{eq:
estimating equation} 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
\begin{align*}
&\left.\frac{\partial}{\partial\boldsymbol{\gamma}}\mathrm{E}[m'_1(D\beta_0 + \boldsymbol{W}^{\top}\boldsymbol{\gamma},\boldsymbol{Y})(D-\boldsymbol{W}^{\top}\boldsymbol{\mu}_0)]\right|_{\boldsymbol{\gamma}=\boldsymbol{\gamma}_0} = \mathbf 0_{p-1}\quad\text{and}\\
&\left.\frac{\partial}{\partial\boldsymbol{\mu}}\mathrm{E}[m'_1(D\beta_0 + \boldsymbol{W}^{\top}\boldsymbol{\gamma}_0,\boldsymbol{Y})(D-\boldsymbol{W}^{\top}\boldsymbol{\mu})]\right|_{\boldsymbol{\mu}=\boldsymbol{\mu}_0} = \mathbf 0_{p-1},
\end{align*}
which follows from \eqref{eq: definition of mu0} and \eqref{eq:EstimandIntro},
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:
\begin{alg}[\textbf{Three-Step Debiasing}]\label{alg:ThreeStepDebiasing}
Given rules for choosing penalty levels
$\lambda_1,\lambda_2\in(0,\infty)$,\footnote{Here, $\lambda_1$ can be chosen via the BCV method, and
$\lambda_2$ can be chosen either via the BCV method or via the SNMD theory for weighted LASSO, as
discussed in \cite{belloni2016post}.} follow the steps below to obtain a debiased estimator
$\widehat\beta$ of $\beta_0$:
\begin{itemize}
\item[\emph{Step 1 (Initiate):}]
\begin{inparaenum}[a.]
\item Define the (preliminary) estimator $\widetilde\boldsymbol{\theta} =
(\widetilde\beta,\widetilde\boldsymbol{\gamma}^{\top})^{\top}$~of $\boldsymbol{\theta}_0 =
(\beta_0,\boldsymbol{\gamma}_0^{\top})^{\top}$ by
$$
\widetilde\boldsymbol{\theta} \in \operatornamewithlimits{argmin}\limits_{\boldsymbol{\theta}\in\Theta}\left\{\mathbb{E}_{n}[m(\boldsymbol{X}_i^{\top}\boldsymbol{\theta},\boldsymbol{Y}_i)] + \lambda_1\|\boldsymbol{\theta}\|_1\right\}.
$$
\item (Optional): Define $\widetilde{T}_1:=\mathrm{supp}(\widetilde\boldsymbol{\theta})$ and recast $\widetilde{\boldsymbol{\theta}}$ as the refitted estimator of $\boldsymbol{\theta}_0$:
$$
\widetilde\boldsymbol{\theta} \in \operatornamewithlimits{argmin}\limits_{\mathclap{\substack{\boldsymbol{\theta}\in\Theta,\\\mathrm{supp}(\boldsymbol{\theta})\subseteq\widetilde{T}_1}}}\mathbb{E}_{n}[m(\boldsymbol{X}_i^{\top}\boldsymbol{\theta},\boldsymbol{Y}_i)].
$$
\end{inparaenum}
\item[\emph{Step 2 (Orthogonalize):}]
\begin{inparaenum}[a.]
\item Based on $\widetilde{\boldsymbol{\theta}}$ from Step 1, define an estimator $\widetilde\boldsymbol{\mu}$ of $\boldsymbol{\mu}_0$ by
\begin{equation}\label{eq: estimator mu}
\widetilde\boldsymbol{\mu} \in \operatornamewithlimits{argmin}\limits_{\boldsymbol{\mu}\in\mathbb{R}^{p-1}}\left\{\mathbb{E}_{n}\big[ m''_{11}(\boldsymbol{X}_i^{\top}\widetilde\boldsymbol{\theta},\boldsymbol{Y}_i)(D_i - \boldsymbol{W}_i^{\top}\boldsymbol{\mu})^2\big] + \lambda_2\|\boldsymbol{\mu}\|_1\right\}.
\end{equation}
\item (Optional): Define $\widetilde{T}_2:=\mathrm{supp}(\widetilde\boldsymbol{\mu})$ and recast $\widetilde{\boldsymbol{\mu}}$ as the refitted estimator of $\boldsymbol{\mu}_0$:
\begin{equation}\label{eq: estimator mu post}
\widetilde\boldsymbol{\mu} \in \operatornamewithlimits{argmin}\limits_{\mathclap{\substack{\boldsymbol{\mu}\in\mathbb{R}^{p-1},\\\mathrm{supp}(\boldsymbol{\mu})\subseteq\widetilde{T}_2}}}\mathbb{E}_{n}\big[ m''_{11}(\boldsymbol{X}_i^{\top}\widetilde\boldsymbol{\theta},\boldsymbol{Y}_i)(D_i - \boldsymbol{W}_i^{\top}\boldsymbol{\mu})^2\big].
\end{equation}
\end{inparaenum}
\item[\emph{Step 3 (Update):}] Define the (debiased) estimator $\widehat\beta$ of $\beta_0$ as the
one-step update of $\widetilde\beta$:
\begin{equation}\label{eq: debiased estimator betahat}
\widehat\beta := \widetilde\beta - \frac{\mathbb{E}_{n}\big[m'_1(\boldsymbol{X}_i^{\top}\widetilde\boldsymbol{\theta}, \boldsymbol{Y}_i)(D_i - \boldsymbol{W}_i^{\top}\widetilde\boldsymbol{\mu})\big]}{\mathbb{E}_{n}\big[m''_{11}(\boldsymbol{X}_i^{\top}\widetilde\boldsymbol{\theta}, \boldsymbol{Y}_i)(D_i - \boldsymbol{W}_i^{\top}\widetilde\boldsymbol{\mu})D_i\big]}.\footnote{In Examples \ref{exa:Expectile} and \ref{exa:PanelCensoredRegressionAndTrimming}, the second derivative $m_{11}''(\boldsymbol{X}_i^{\top}\widetilde\boldsymbol{\theta},\boldsymbol{Y}_i)$ may not exist for some observations $i\in[n]$, in which case the estimator $\widehat\beta$ may not be well-defined. To make it well-defined, we replace $m_{11}''(\boldsymbol{X}_i^{\top}\widetilde\boldsymbol{\theta},\boldsymbol{Y}_i)$ in \eqref{eq: estimator mu}, \eqref{eq: estimator mu post} and \eqref{eq: debiased estimator betahat} for such observations by zero.}
\end{equation}
\end{itemize}
\end{alg}
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.
\begin{assumption}[\textbf{Identifiability}]\label{as: identifiability inference} There exists a constant $c_I\in(0,\infty)$ such that we have
$\mathrm{E}[|m_1'(\boldsymbol{X}^{\top}\boldsymbol{\theta}_0,\boldsymbol{Y})(D-\boldsymbol{W}^{\top}\boldsymbol{\mu}_0)|^2]\geqslant c_I$.
\end{assumption}
\begin{assumption}[\textbf{Integrability}]\label{as: integrability inference}
There are constants $C_M\in(0,\infty)$ and $\widetilde r\in(4,\infty)$ such that
$\max_{1\leqslant j\leqslant p}(\mathrm{E}[|X_{j}|^{\widetilde r}])^{1/\widetilde
r}\leqslant C_M$, $(\mathrm{E}[|D - \boldsymbol{W}^{\top}\boldsymbol{\mu}_0|^{\widetilde r}])^{1/\widetilde
r}\leqslant C_M$, and $(\mathrm{E}[|m_1'(\boldsymbol{X}^{\top}\boldsymbol{\theta}_0,\boldsymbol{Y})|^{\widetilde
r}])^{1/\widetilde r}\leqslant C_M$.
\end{assumption}
Assumption \ref{as: identifiability inference} 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{sec:VerificationLinearModelSquareLoss} for
details.\footnote{More generally, Assumption \ref{as: identifiability inference}
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{as: integrability inference} 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.
\begin{assumption}[\textbf{Smoothness}]\label{as: smoothness inference} There
are constants $C_m\in(0,\infty)$, $J\in\mathbb N$, and a possibly $\boldsymbol{y}$-dependent
partition $-\infty=t_{\boldsymbol{y},0} < t_{\boldsymbol{y},1} \leqslant \dots \leqslant t_{\boldsymbol{y},J-1} <
t_{\boldsymbol{y},J} = \infty$ of $\mathbb R$ such that for all $\boldsymbol{y}\in\mathcal Y$, the
function $m(\cdot,\boldsymbol{y})$ is continuously differentiable on $\mathbb{R}$ and three-times
differentiable on each $(t_{\boldsymbol{y},j-1},t_{\boldsymbol{y},j})$, $j\in[J]$, with second and
third derivatives satisfying $|m''_{11}(t,\boldsymbol{y})|\leqslant C_m$ and
$|m'''_{111}(t,\boldsymbol{y})|\leqslant C_m$. In addition,
$m_{11}''(\boldsymbol{X}^{\top}\boldsymbol{\theta}_0,\boldsymbol{Y})$ exists almost surely.
\end{assumption}
This assumption strengthens Assumption \ref{as: post estimator smoothness} from Section
\ref{sec:Bootstrapping-the-Penalty} (which is why we reuse the symbol $C_m$ for the constant). Note
that Assumption \ref{as: smoothness inference} does not hold for the trimmed LAD loss function in
Example \ref{exa:PanelCensoredRegressionAndTrimming}, which means that our inference approach does not apply for this loss function. In addition, Assumption \ref{as: smoothness inference} 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 \cite{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{sec: verification}, we verify Assumption \ref{as: smoothness inference} for all other
examples from Section \ref{sec:Examples} including Example
\ref{exa:PanelCensoredRegressionAndTrimming} 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.
\begin{assumption}[\textbf{Density}]\label{as: conditional density}
Provided $J\geqslant 2$, there is a constant $C_f\in(0,\infty)$ and a non-random
sequence $\overline\Delta_n$ in $(0,\infty)$ such that $\mathrm{P}(\boldsymbol{X}^{\top}\boldsymbol{\theta}_0 -
\overline\Delta_n \leqslant t_{\boldsymbol{Y},j}\leqslant \boldsymbol{X}^{\top}\boldsymbol{\theta}_0 +
\overline\Delta_n) \leqslant C_f\overline\Delta_n$ for all
$j\in[J-1]$.
\end{assumption}
Assumption \ref{as: conditional density} holds trivially in Examples \ref{exa:Logit} and
\ref{exa:LogConcaveOrderedResponse}, as for those examples Assumption \ref{as: smoothness inference}
holds with $J=1$. When combined with the requirement that $\overline\Delta_n\to0$ (sufficiently
fast), Assumption \ref{as: conditional density} does impose quite a bit of structure in Examples
\ref{exa:Expectile} and \ref{exa:PanelCensoredRegressionAndTrimming}, however. In Example
\ref{exa:Expectile}, this assumption is satisfied if the conditional distribution of $Y$ given $\boldsymbol{X}$
is absolutely continuous with bounded PDF. In Example \ref{exa:PanelCensoredRegressionAndTrimming}
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{sec: verification} 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$.}
\begin{assumption}[\textbf{Convergence Rates}]\label{as: convergence rates inference}
There is a non-random sequence $a_n$ in $(0,\infty)$ such that
$a_n\to 0$ and $\|\widetilde\boldsymbol{\theta} - \boldsymbol{\theta}_0\|_1 + \|\widetilde\boldsymbol{\mu} -
\boldsymbol{\mu}_0\|_1 \lesssim_{\mathrm{P}} a_n$.
\end{assumption}
Assumption \ref{as: convergence rates inference} is a high-level assumption placed on the estimators
from Steps 1 and 2 in Algorithm \ref{alg:ThreeStepDebiasing}. 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{thm: convergence rates post estimator}. For the estimation error of
$\widetilde\boldsymbol{\mu}$, however, we cannot use Theorem \ref{thm: convergence rates post estimator}, as
this estimator does not fit into our framework because of the presence of estimated weights in the
optimization problems \eqref{eq: estimator mu} and \eqref{eq: estimator mu post}. 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.~\cite{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$.
\begin{thm}[\textbf{Asymptotic Distribution}]\label{thm: asymptotic distribution} Let Assumptions
\ref{assu:ParameterSpace}--\ref{assu:Approximate-Sparsity} and \ref{as: identifiability
inference}--\ref{as: convergence rates inference} hold, and suppose that $\sqrt n a_n^2\to0$,
$a_n(n^{1/r}B_n + \sqrt{\ln(pn)})\to 0$, and $B_n^2\ln(pn) = o(n^{1-4/(r\wedge\widetilde r)})$.
If $J\geqslant 2$, suppose also that $(1+\sqrt n B_na_n)(\overline\Delta_n^{1/2} +
(B_na_n/\overline\Delta_n)^{r/2})\to0$. Then
\begin{equation}\label{eq: asy normality}
\frac{\sqrt n(\widehat\beta - \beta_0)}{\sigma_{0}}\overset{D}\to \mathrm{N}(0,1),\quad\text{where}\quad \sigma_{0}^2:=\frac{\mathrm{E}\big[(m'_1(\boldsymbol{X}^{\top}\boldsymbol{\theta}_0,\boldsymbol{Y})(D-\boldsymbol{W}^{\top}\boldsymbol{\mu}_0))^2\big]}{\big(\mathrm{E}[m''_{11}(\boldsymbol{X}^{\top}\boldsymbol{\theta}_0,\boldsymbol{Y})(D-\boldsymbol{W}^{\top}\boldsymbol{\mu}_0)D]\big)^2}.
\end{equation}
\end{thm}
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
\begin{equation}\label{eq: sigma hat 1}
\widehat{\sigma}^2:=\frac{\mathbb{E}_{n}\big[(m'_1(D_i\widetilde\beta + \boldsymbol{W}_i^{\top}\widetilde\boldsymbol{\gamma}, \boldsymbol{Y}_i)(D_i - \boldsymbol{W}_i^{\top}\widetilde\boldsymbol{\mu}))^2\big]}{\big(\mathbb{E}_{n}\big[m''_{11}(D_i\widetilde\beta + \boldsymbol{W}_i^{\top}\widetilde\boldsymbol{\gamma}, \boldsymbol{Y}_i)(D_i - \boldsymbol{W}_i^{\top}\widetilde\boldsymbol{\mu})D_i\big]\big)^2}
\end{equation}
with the estimators $\widetilde\beta,\widetilde\boldsymbol{\gamma}$ and $\widetilde\boldsymbol{\mu}$ stemming from Steps 1
and 2 of Algorithm \ref{alg:ThreeStepDebiasing} (possibly with refitting) using BCV as the
penalty rule in both steps. Alternatively, one can incorporate Step 3 of the same algorithm and use
\begin{equation}\label{eq: sigma hat 2}
\widehat{\sigma}^2:=\frac{\mathbb{E}_{n}\big[(m'_1(D_i\widehat\beta + \boldsymbol{W}_i^{\top}\widetilde\boldsymbol{\gamma}, \boldsymbol{Y}_i)(D_i - \boldsymbol{W}_i^{\top}\widetilde\boldsymbol{\mu}))^2\big]}{\big(\mathbb{E}_{n}\big[m''_{11}(D_i\widehat\beta + \boldsymbol{W}_i^{\top}\widetilde\boldsymbol{\gamma}, \boldsymbol{Y}_i)(D_i - \boldsymbol{W}_i^{\top}\widetilde\boldsymbol{\mu})D_i\big]\big)^2}.\footnote{For both variance estimators \eqref{eq: sigma hat 1} and \eqref{eq: sigma hat 2}, in case the second derivative $m_{11}''(\boldsymbol{X}_i^{\top}\boldsymbol{\theta},\boldsymbol{Y}_i)$ fails to exist at $\boldsymbol{\theta}=(\widetilde\beta,\widetilde\boldsymbol{\gamma}^\top)^\top$ or $\boldsymbol{\theta}=(\widehat\beta,\widetilde\boldsymbol{\gamma}^\top)^\top$ and some $i\in[n]$, we replace those second derivatives by zero.}
\end{equation}
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 \eqref{eq: sigma hat
1} or by \eqref{eq: sigma hat 2}, and $z_{\alpha/2}$ is the $(1-\alpha/2)$-quantile of the standard
normal distribution.
\begin{rem}[{\textbf{Relation to Literature}}]
As discussed in the beginning of this section, our approach to inference in this
section closely follows the developments in the literature. In particular, our
estimator $\widehat\beta$ is essentially the same as that proposed in
\cite{van_de_geer_high-dimensional_2014}, the only difference being that we
allow refitting in the optional parts of Steps 1 and 2 of Algorithm
\ref{alg:ThreeStepDebiasing}. As we will see in the next section, this refitting
can substantially improve inference, in terms of size control, even in
approximately sparse models. More importantly, however, is that Theorem
\ref{thm: asymptotic distribution} is different from the corresponding theorem
in \cite{van_de_geer_high-dimensional_2014}, as we tune the assumptions of our
theorem toward the examples from Section \ref{sec:Examples}. Specifically, we do
not require the function $m(\cdot,\boldsymbol{y})$ to be strictly convex (see \textit{ibid.},
p.~1179) or for it to be everywhere twice differentiable with a
Lipschitz-continuous second derivative (see \textit{ibid.}, Condition (C1)). No
matter the choice of ``trimmer'' $\Xi$, the trimmed loss \eqref{eq:TrimmedLoss} in
Example \ref{exa:PanelCensoredRegressionAndTrimming} has linear pieces and is
therefore not strictly convex. Moreover, neither the asymmetric LS (Example
\ref{exa:Expectile} with $\tau\neq1/2$) nor trimmed LS loss functions have
Lipschitz-continuous second derivatives. We also provide a detailed verification
of our assumptions for all Section \ref{sec:Examples} examples in Appendix
\ref{sec: verification}. Related approaches to debiasing of high-dimensional
estimators were also proposed in \cite{javanmard_confidence_2013} in a
likelihood framework and in \cite{belloni2016post} for generalized linear
models.\qed
\end{rem}
\section{Simulations}\label{sec: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{sec:Bootstrapping-the-Penalty}. 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.
\subsection{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{exa:Logit}. 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:
\begin{align*}
\mathrm{Pattern\,1}:\quad & \boldsymbol{\theta}_{0}=\left(1,1,0,\dotsc,0\right)^{\top},\tag{Exactly Sparse}\\
\mathrm{Pattern\,2}:\quad & \theta_{0,j}=(1/\sqrt{2})^{j-1}\mathbf{1}\left(j\leqslant5\right),\quad j\in\left[p\right],\tag{Intermediate}\\
\mathrm{Pattern\,3}:\quad & \theta_{0,j}=(1/\sqrt{2})^{j-1},\quad j\in\left[p\right].\tag{Approximately Sparse}
\end{align*}
The \emph{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
\begin{equation}
\mathrm{var}\left(\boldsymbol{X}^{\top}\boldsymbol{\theta}_{0}\right)=2\left(1+\rho\right)\in\left\{ 2,2.4,2.8,3.2,3.8\right\} .\label{eq:SimulationsSignal}
\end{equation}
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 \citet[Section 5.2]{friedman_regularization_2010} and
\citet[Section 5]{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
\emph{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 \emph{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{assu:Approximate-Sparsity}, 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 \emph{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.
\begin{rem}[\textbf{Sparsity of Debiasing Coefficient Vector}]
One may wonder whether the above patterns for the structural coefficients
$\boldsymbol{\theta}_{0}=(\beta_{0},\boldsymbol{\gamma}_{0}^{\top})^{\top}$ agree or conflict with sparsity of the
non-primitive debiasing coefficient vector $\boldsymbol{\mu}_{0}$ in any sense of the word. While a thorough
investigation of this question is beyond the scope of this paper, we can provide some insights for
our concrete DGPs. Specifically, in Appendix \ref{sec:Sparsity-of-Debiasing-Vector}, we show
that---in our collection of DGPs---the number of non-zeros in $\boldsymbol{\mu}_{0}$ is bounded by the number of
non-zeros in $\boldsymbol{\gamma}_{0}$, $\|\boldsymbol{\mu}_{0}\|_{0}\leqslant\|\boldsymbol{\gamma}_{0}\|_{0}$, i.e.~the number of
relevant controls. Hence, when $\boldsymbol{\gamma}_0$ is (exactly) sparse, so is $\boldsymbol{\mu}_0$. Moreover, we show
via simulation that when $\boldsymbol{\gamma}_{0}$ is only approximately sparse, the sorted absolute values of
the elements of $\boldsymbol{\mu}_{0}$ are rapidly decaying and approaching zero, cf.~Figure
\ref{fig:SortedAbsoluteDebiasingVector}. Such a decay is in line with the notion of approximate
sparsity.\qed
\end{rem}
\subsection{Estimation and Implementation}
We consider the following four estimators arising from $\ell_1$-ME
\eqref{eq:ell1PenalizedMEstimationIntro} and post-$\ell_1$-ME \eqref{eq:AllPostEstimators} based on
either the CV or BCV penalty levels in \eqref{eq:LambdaCVDefn} and \eqref{eq:BootstrapPenaltyLevel
CV}, respectively:
\begin{itemize}
\item $\ell_1$-ME based on bootstrapping after cross-validation (``BCV''),
\item post-$\ell_1$-ME based on bootstrapping after cross-validation (``post-BCV''),
\item $\ell_1$-ME based on cross-validation (``CV''), and
\item post-$\ell_1$-ME based on cross-validation (``post-CV'').
\end{itemize}
When discussing normal approximations based on three-step debiasing (Algorithm
\ref{alg:ThreeStepDebiasing}), 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 \citet[p.~2380]{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{sec:Bootstrapping-the-Penalty}. The alternative probability tolerance rule
$\alpha_{n}=10/n$ leads to qualitatively identical conclusions, cf.~Appendix
\ref{sec:Additional-Simulation-Results}. 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{alg:ThreeStepDebiasing} 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
\texttt{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 \texttt{R} with cross-validation done using
\texttt{glmnet::cv.glmnet}, and refitting done using
\texttt{stats::glm}.\footnote{We use \texttt{R} version 4.2.2 and \texttt{glmnet
}version 4.1-6.} When constructing the candidate penalty set $\Lambda_{n}$, we
use the \texttt{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{assu:CandidatePenalties} (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{assu:CandidatePenalties}. 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 \texttt{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.
\subsection{Results}
\subsubsection{Non-Existence and Treatment of Missing Values}
While the $\ell_{1}$-penalized probit estimators BCV and CV always exist (cf.~Section
\ref{sec:ExistenceSparsityAndUniqueness}), 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{sec:ExistenceSparsityAndUniqueness} and, in
particular, Remark \ref{rem:UnpenalizedCoefficients} 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.
\subsubsection{Estimation Error}
Figure \ref{fig:MeanEll2ErrorEstimatorComparisonBCCHTuning} 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.
\begin{figure}
\caption{Mean $\ell_{2}$ Estimation Error by Method with $c_{0}=1.1$ and
$\alpha_{n}=.1/\ln\left(p\lor
n\right)$\label{fig:MeanEll2ErrorEstimatorComparisonBCCHTuning}}
\centering{}\includegraphics[viewport=5bp 5bp 463bp
416bp,clip,width=0.7\textwidth]{Fig_probit_mean_ell2_err_BCV_PostBCV_CV_R2000_B1000_K3_c1dot1_alphabcch_interceptTRUE_standardizeTRUE}
\end{figure}
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 \emph{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{cor: convergence rate bootstrap after cv} and \ref{thm: convergence rates post estimator}.
We next investigate the impact of the choice of score markup $c_{0}$. Figures
\ref{fig:MeanEll2ErrorBCVVaryingc0BCCHrule} and
\ref{fig:MeanEll2ErrorPostBCVVaryingc0BCCHrule} 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.
\begin{figure}[p]
\caption{Mean $\ell_{2}$ BCV Estimation Error by Score Markup with
$\alpha_{n}=.1/\ln\left(p\lor
n\right)$\label{fig:MeanEll2ErrorBCVVaryingc0BCCHrule}}
\centering{}\includegraphics[viewport=5bp 5bp 463bp
416bp,clip,width=0.7\textwidth]{Fig_probit_mean_ell2_BCV_err_for_different_c0_R2000_B1000_K3_alphabcch_interceptTRUE_standardizeTRUE}
\end{figure}
\begin{figure}
\caption{Mean $\ell_{2}$ Post-BCV Estimation Error by Score Markup with
$\alpha_{n}=.1/\ln\left(p\lor
n\right)$\label{fig:MeanEll2ErrorPostBCVVaryingc0BCCHrule}}
\centering{}\includegraphics[viewport=5bp 5bp 463bp
416bp,clip,width=0.7\textwidth]{Fig_probit_mean_ell2_Post-BCV_err_for_different_c0_R2000_B1000_K3_alphabcch_interceptTRUE_standardizeTRUE}
\end{figure}
Figure \ref{fig:MeanEll2ErrorBCVVaryingc0BCCHrule} 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{fig:MeanEll2ErrorPostBCVVaryingc0BCCHrule}), 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.
\subsubsection{Normal Approximation}
We next assess the normal approximations resulting from three-step debiasing
(Algorithm \ref{alg:ThreeStepDebiasing}) using either BCV, post-BCV or CV.
Instead of looking at the \emph{standardized} estimate
$\sqrt{n}(\widehat{\beta}-\beta_{0})/\sigma_{0}$ for the true asymptotic
variance $\sigma_0^2$ given in \eqref{eq: asy normality},
we form an estimate $\widehat{\sigma}^{2}$ and consider the \emph{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
\begin{align*}
\mathrm{E}\left[ m_{11}''\left(\boldsymbol{X}^{\top}\boldsymbol{\theta}_{0},\boldsymbol{Y}\right) \mid \boldsymbol{X}^\top\boldsymbol{\theta}_0 = t \right]
= \frac{f\left( t \right)^2 }{F\left( t \right)\left( 1 - F \left( t \right) \right)}
=: \omega_{F}\left(t\right),
\end{align*}
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{sec:Application}, 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{alg:ThreeStepDebiasing} given
different rules for choosing the penalties $\lambda_1$ and
$\lambda_2$.\footnote{Alternatively, one can use the
``sandwich'' estimators \eqref{eq: sigma hat 1} and \eqref{eq: sigma hat 2}.
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{fig:NormalApproxBCVPostBCVCVVaryingSampleSizeBCCHtuningZeroCorr}
shows the (kernel) densities of the studentized estimates using
benchmark tuning and $\rho=0$.\footnote{All kernel densities are created using
the \texttt{R} package \texttt{ggplot2} with \texttt{geom\_density}. In
expectation of an approximately normal distribution, we use a Gaussian kernel
and the \citet[Equation (3.31)]{silverman_density_1986} rule-of-thumb bandwidth
(both \texttt{geom\_density }defaults).}
\begin{figure}
\caption{{\small{}Densities of Studentized Estimates by $n(=p)$ with $\rho=0$,
$c_{0}=1.1$ and $\alpha_{n}=.1/\ln\left(p\lor
n\right)$.}\label{fig:NormalApproxBCVPostBCVCVVaryingSampleSizeBCCHtuningZeroCorr}}
\centering{}\includegraphics[viewport=5bp 5bp 463bp
416bp,clip,width=0.68\textwidth]{Fig_probit_normal_approx_by_n_each_method_and_pattern_R2000_B1000_K3_c1dot1_alphabcch_rhodot0_interceptTRUE_standardizeTRUE}
\end{figure}
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{fig:NormalApproxBCVPostBCVCVVaryingSampleSizeBCCHtuningZeroCorr} depicts
the normal approximations for the worst-correlation case $\rho=0$, in Figure
\ref{fig:NormalApproxPostBCVCVBCCHtuningVaryingCorr} we display the normal
approximations as a function of $\rho$.
\begin{figure}[!htb]
\caption{{\footnotesize{}Densities of Studentized Post-BCV and CV Estimates for
Different $\rho$ with $n(=p)=400,$ Approximately Sparse Coefficient Pattern, $c_{0}=1.1$ and
$\alpha_{n}=.1/\ln\left(p\lor
n\right)$}\label{fig:NormalApproxPostBCVCVBCCHtuningVaryingCorr}}
\centering{}\includegraphics[viewport=5bp 5bp 463bp
416bp,clip,width=0.68\textwidth]{Fig_probit_normal_approx_Post-BCV_vs_CV_for_different_rho_R2000_B1000_K3_dense_np400_c1dot1_alphabcch_interceptTRUE_standardizeTRUE}
\end{figure}
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{fig:MeanEll2ErrorPostBCVVaryingc0BCCHrule}, the exact markup value
appears to make little difference. Figures are available upon request.}
\section{Revisiting Racial Differences in Police Use of Force}\label{sec:Application}
In this section we revisit the empirical setting in
\citet{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{exa:Logit}.\footnote{See \citet{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
\begin{equation}
\mathrm{P}\left(\mathrm{Force}=1\mid\mathbf{Race},\boldsymbol{W}\right)=F\left(\mathbf{Race}^{\top}\boldsymbol{\alpha}_{0}+\boldsymbol{W}^{\top}\boldsymbol{\gamma}_{0}\right),\label{eq:FryerModel}
\end{equation}
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 \citet[Table
2.B]{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 \citet[Table
2.B]{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
\begin{inparaenum}[(i)]
\item an alternative binary response model and
\item large(r) sets of candidate regressors, in combination with
$\ell_{1}$-penalization.
\end{inparaenum}
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 \eqref{eq:FryerModel} 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 \citet{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 \citet[Table
2.B]{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{alg:ThreeStepDebiasing}) 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 \citet[Table 2.B, Row
l]{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 \emph{total} number of police-civilian encounters in the PPCS dataset
covering the six surveys 1996, 1999, 2002, 2005, 2008 and 2011. The number of \emph{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{tab:RevisitingFryerTstats} displays the $t$-values associated
with testing the null hypothesis using either unpenalized or
$\ell_{1}$-penalized methods.
\begin{table}
\caption{$t$-Values for Testing a Zero Coefficient
$\beta_{\mathrm{Black}}$ on $\mathrm{Black}$\label{tab:RevisitingFryerTstats}}
\centering{}
\begin{tabular}{lrrrrr}
\hhline{======}
& \multicolumn{2}{c}{Unpenalized (ML)} & &
\multicolumn{2}{c}{Post-BCV}\tabularnewline \cline{2-3} \cline{3-3} \cline{5-6}
\cline{6-6} Controls \textbackslash{} Loss & Logit & Probit & & Logit &
Probit\tabularnewline
\hline
Basic Controls & 8.8 & 8.7 & & 10.5 & 9.6\tabularnewline + Interactions & n.a.
& n.a. & & 20.7 & 18.9\tabularnewline
\hline
\end{tabular}
\end{table}
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{tab:RevisitingFryerTstats}
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
\citet{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{eq:FryerModel}), the APE can be
expressed as the average probability difference
\begin{align}
\mathrm{APE}_{\mathrm{Black}} & :=\mathrm{E}\big[\mathrm{P}\left(\mathrm{Force}=1\mid\mathrm{Black}=1,\mathrm{Hisp}=0,\mathrm{Other}=0,\boldsymbol{W}\right)\nonumber \\
& \quad\quad-\mathrm{P}\left(\mathrm{Force}=1\mid\mathrm{Black}=0,\mathrm{Hisp}=0,\mathrm{Other}=0,\boldsymbol{W}\right)\big]\nonumber \\
& =\mathrm{E}\left[F\left(\beta_{\mathrm{Black}}+\boldsymbol{W}^{\top}\boldsymbol{\gamma}_{0}\right)-F\left(\boldsymbol{W}^{\top}\boldsymbol{\gamma}_{0}\right)\right],\label{eq:APE}
\end{align}
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{tab:RevisitingFryerAPEs} reports the APE
estimates (in percentage points) corresponding to these procedures including
either basic controls or basic controls with interactions.
\begin{table}[htb]
\caption{Estimates of the Average Partial Effect of $\text{Black}$ (in Percentage
Points)\label{tab:RevisitingFryerAPEs}}
\centering{}
\begin{tabular}{lrrrrr}
\hhline{======}
& \multicolumn{2}{c}{Unpenalized (ML)}
& & \multicolumn{2}{c}{Post-BCV}\tabularnewline \cline{2-3} \cline{3-3}
\cline{5-6} \cline{6-6} Controls \textbackslash{} Loss & Logit & Probit & &
Logit & Probit\tabularnewline
\hline
Basic Controls & 1.1 & 1.1 & & 1.4 & 1.3\tabularnewline + Interactions & n.a. &
n.a. & & 3.2 & 2.8\tabularnewline
\hline
\end{tabular}
\end{table}
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{eq:APE}) 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{tab:Timings} we
report the computing time used by the above-mentioned estimation routines.
\begin{table}[htb]
\begin{threeparttable}
\caption{Estimation Routine Timings (in Seconds) \label{tab:Timings}}
\centering{}
\begin{tabular}{lrrrrr}
\hhline{======}
& \multicolumn{2}{c}{Unpenalized (ML)} & &
\multicolumn{2}{c}{Post-BCV}\tabularnewline \cline{2-3} \cline{3-3}
\cline{5-6} \cline{6-6} Controls \textbackslash{} Loss & Logit & Probit &
& Logit & Probit\tabularnewline
\hline
Basic Controls & 1.5 & 1.5 & & 21 & 40\tabularnewline + Interactions
& $\infty$ & $\infty$ & & 210 & 199\tabularnewline
\hline
\end{tabular}
\begin{tablenotes}
\footnotesize
\item\emph{Notes:} All timings were carried out on an Intel Core
i7-8700 3.20GHz CPU. When using \texttt{cv.glmnet}, we use the parallel
computing option with all 12 virtual cores available.
\end{tablenotes}
\end{threeparttable}
\end{table}
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 \texttt{glmnet} does \emph{not} check for optimizer
divergence \citep[cf.][p.~9]{friedman_regularization_2010}. We therefore opted for
\texttt{stats::glm} for ML estimation and refitting.}
Of course, the Table \ref{tab:Timings} runtimes are only single observations arising from our
particular \texttt{R} implementation of our procedures, using a specific dataset, and our specific
computing environment. As such, they need not translate to other settings.
\bibliographystyle{ecta}
\bibliography{My_Library}
\newpage{}