# Debiased Inference on Heterogeneous Quantile Treatment Effects with
Regression Rank-Scores

Alexander Giessing Department of Statistics, University of Washington. E-mail:
<EMAIL_ADDRESS>Jingshen Wang Division of Biostatistics, University of
California, Berkeley. E-mail<EMAIL_ADDRESS>
###### Abstract

Understanding treatment effect heterogeneity is vital to many scientific
fields because often the same treatment affects different individuals
differently. Quantile regression provides a natural framework for modeling
such heterogeneity. We propose a new method for inference on heterogeneous
quantile treatment effects in the presence of high-dimensional covariates. Our
estimator combines an $\ell_{1}$-penalized regression adjustment with a
quantile-specific bias correction scheme based on rank scores. We study the
theoretical properties of this estimator, including weak convergence of the
heterogeneous quantile treatment effect process to a Gaussian process. We
illustrate the finite-sample performance of our approach through simulations
and an empirical example, dealing with the differential effect of statin usage
for lowering low-density lipoprotein cholesterol levels for the Alzheimer’s
disease patients who participated in the UK Biobank study.

Keywords: Quantile Regression; Debiased Inference; High-dimensional Data;
Semiparametric Efficiency; Causal Inference.

##  1 Introduction

###  1.1 Motivation

Understanding treatment effect heterogeneity in observational studies is vital
to many scientific fields because often the same treatment affects different
individuals differently. For instance, in modern drug development, it is
important to test for the existence (or the lack) of treatment effect
heterogeneity and to identify subpopulations for which a treatment is most
beneficial (or harmful) (Lipkovich et al.,, 2017; Ma and Huang,, 2017).
Similarly, in precision medicine, it is essential to be able to generalize
causal effect estimates from a small experimental sample to a target
population (Kern et al.,, 2016; Coppock et al.,, 2018). Since quantile
regression (Koenker,, 2005) models the effect of covariates on the conditional
distribution of the response variable, it provides a natural framework for
studying treatment heterogeneity. In this paper, we propose a new procedure
for inference on the heterogeneous quantile treatment effects (HQTE) curve in
the presence of high-dimensional covariates. The HQTE curve is defined as the
difference between the quantiles of the conditional distributions of treatment
and control group:

$\displaystyle\alpha(\tau;z):=Q_{1}(\tau;z)-Q_{0}(\tau;z),$ (1)

where $Q_{1}(\tau;z)$ ($Q_{0}(\tau;z)$) is the conditional quantile curve of
the potential outcome of the treated group (the control group) evaluated at a
quantile level $\tau\in(0,1)$ and covariate $z\in\mathbb{R}^{p}$.

The HQTE is a versatile quantity to model treatment effect heterogeneity
because the HQTE curve provides information about the treatment effect at
every quantile level. It therefore affords a more nuanced analysis of
treatment effects than the average treatment effect, which only provides
information about the mean of the outcome variable. For example, in the
canonical study of infant birth weights, literature documents that maternal
hypertension is a risk factor for lower birth weight and that the effect of
hypertension on birth weight is greater on the left tail of the birth weight
distribution (Bowers et al.,, 2011; Mhanna et al.,, 2015). In this case, a
statistical procedure aimed at detecting treatment effects at the lower
quantiles of the birth weight distribution is more useful than a statistical
procedure designed to estimate the average treatment effect. As another
example, when the outcome variable has a skewed distribution (such as survival
times), the HQTE curve is a more informative quantity than the classical
average treatment effect. Wang et al., (2018) show that mean-optimal treatment
regimes that maximize the average treatment effect can be harmful to
individuals who differ significantly from the typical individual in the
sample. In such a case, a more adaptive quantile-optimal treatment regime
based on the HQTE is preferred.

###  1.2 Our contribution

The primary contribution of this article is the novel _rank-score debiased
estimator_ of the heterogeneous quantile treatment effects (HQTE) curve and a
comprehensive study of its theoretical properties. We break down our
contribution as follows:

On the methodological side, we show how to use inverse-density weighted
regression rank-scores to debias estimates of the conditional quantile
function when these estimates are obtained from solving an
$\ell_{1}$-penalized quantile regression problem. We rationalize this idea in
two different ways: a bias-variance trade-off and an approximate Neyman
orthogonalization procedure (Section 3).

On the theoretical side, our main result is the weak convergence of the rank-
score debiased HQTE curve to a Gaussian process in
$\ell^{\infty}(\mathcal{T})$. The large sample properties of this process are
needed whenever one would like to conduct simultaneous inference on the HQTE
curve on several (or a continuum) of quantile levels
$\mathcal{T}\subset(0,1)$. We propose two uniformly consistent estimators for
the covariance functions of the Gaussian limit process. Moreover, for fixed
dimensions, we prove that the rank-score debiased estimator is semiparametric
efficient (Section 4).

On the practical side, we propose a systematic way of selecting the tuning
parameters in the proposed estimation procedure. Our procedure is similar to
optimization problems adopted for covariate balancing in causal inference
(Zubizarreta,, 2015; Wang and Zubizarreta,, 2017; Athey et al.,, 2018). While
conventional covariate balancing procedures are rather sensitive to the choice
of tuning parameters, our systematic procedure makes use of the dual
formulation of the rank-score debiasing program and is fully automatic
(Section 5). We illustrate the finite-sample performance of our approach
through Monte Carlo experiments (Section 6) and an empirical example, dealing
with the differential effect of statin usage on lowering the low-density
lipoprotein cholesterol (LDL) levels for the Alzheimer’s disease patients
(Section 7).

To analyze the theoretical properties of the quantile rank-score debiasing
problem, we develop new technical tools that complement existing results on
the consistency of $\ell_{1}$-penalized quantile regression (Wang and He,,
2021; Belloni et al., 2019b, ; Belloni and Chernozhukov,, 2011) and the weak
convergence of quantile regression processes in growing dimension (Belloni et
al., 2019a, ; Chao et al.,, 2017). Two new results are particularly
interesting: the dual formulation of the rank-score debiasing program and the
Bahadur-type representation for the rank-score debiased estimator (Sections
C–I of the Supplementary Materials).

###  1.3 Prior and related work

Treatment effect heterogeneity is of significant interest in causal inference
and is analyzed from many different angles. Imai and Ratkovic, (2013)
formulate the estimation of heterogeneous mean treatment effects as a variable
selection problem. Angrist, (2004) studies mean treatment effect heterogeneity
through instrumental variables. In recent publications Semenova and
Chernozhukov, (2021), Künzel et al., (2019) and Nie and Wager, (2019) propose
several new meta-learners to estimate conditional average treatment effects.
Firpo, (2007); Frölich and Melly, (2013); Cattaneo, (2010) study (marginal)
quantile treatment effects through modeling inverse propensity scores.
Chernozhukov and Hansen, (2005) and Abadie et al., (2002) show how
instrumental variables can be helpful in identifying conditional quantile
treatment effects in the presence of unmeasured confounding variables. Our
paper contributes to this thriving field by introducing a novel quantile
estimator to address treatment effect heterogeneity.

Several recent articles study the problem of debiased inference for high-
dimensional quantile regression. Belloni et al., 2019b  propose an efficient
debiased estimator of a single quantile regression coefficient using Neyman
orthogonal scores. Bradic and Kolar, (2017) consider the problem of debiasing
the $\ell_{1}$-penalized estimate of the quantile regression vector when the
response is homoscedastic. Zhao et al., (2019) consider the same problem but
propose different estimator that can deal with heteroscedastic responses.
Allowing for heteroscedastic responses is indeed of great practical importance
since the ability to model heteroscedasticity is a key reason for using
quantile regression in the first place. The important conceptual difference
between these approaches and ours is that we treat the solution of the
$\ell_{1}$-penalized quantile regression problem as a nuisance parameter and
directly debias the scalar estimate of the conditional quantile function
$Q_{d}(\tau;z)$. While at least in principle the estimator by Zhao et al.,
(2019) can also be used to construct debiased estimates of the conditional
quantile function, our approach has the following three advantages: First, we
show via heuristic derivations and simulation studies that our estimator is
more efficient. Moreover, we prove that in fixed dimensions our estimator
achieves the semi-parametric efficiency bound. Second, we derive the limiting
distribution of our estimator, whereas this is not possible using the results
by Zhao et al., (2019). Third, our estimator is computationally efficient as
we do not need to estimate the inverse of a high-dimensional covariance matrix
to debias the quantile regression coefficient vector. We delegate a more
detailed comparison of the approaches by Belloni et al., 2019b  and Zhao et
al., (2019) and ours to Section A of the Supplementary Materials.

##  2 Causal framework and identification

Throughout this paper, $Y\in\mathbb{R}$ denotes the response variable,
$D\in\\{0,1\\}$ a binary treatment variable, and $X\in\mathbb{R}^{p}$ a vector
of covariates. Following the framework of Rubin, (1974), we define the causal
effect of interest in terms of so-called potential outcomes: Potential
outcomes describe counterfactual states of the world, i.e. possible responses
if certain treatments were administered. More formally, we index the outcomes
of the response variable $Y$ by the treatment variable $D$ and write $Y_{D}$
for the potential outcomes of $Y$. With this notation, the potential outcome
$Y_{d}$ corresponds to the response that we would observe if treatment $D=d$
was assigned. The causal quantity of interest in this paper is the
heterogeneous quantile treatment effect (HQTE) curve evaluated at covariates
$z\in\mathbb{R}^{p}$,

$\displaystyle\alpha(\tau;z):=Q_{1}(\tau;z)-Q_{0}(\tau;z),$ (2)

where
$Q_{d}(\tau;z)=\inf\left\\{y\in\mathbb{R}:F_{Y_{d}|X}(y|z)\geq\tau\right\\}$
is the conditional quantile function (CQF) of the potential outcome $Y_{d}\mid
X=z$ at a quantile level $\tau\in(0,1)$ and $F_{Y_{d}|X}$ denotes the
corresponding conditional distribution function.

The key challenge in causal inference is that for each individual we only
observe its potential outcome $Y_{D}$ under one of the two possible treatment
assignments $D\in\\{0,1\\}$ but never under both. In other words, the observed
response variable is given as $Y=DY_{1}+(1-D)Y_{0}$. Since the potential
outcomes $Y_{0}$ and $Y_{1}$ are not observed, a priori, it is unclear how to
estimate $Q_{1}(\tau;z)$ and $Q_{0}(\tau;z)$. To make headway, we introduce
the following condition:

######  Condition 1 (Unconfoundedness).

$(Y_{0},Y_{1})$ is independent of $D$ given $X$, i.e.
$(Y_{0},Y_{1})\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}D\mid X$.

Colloquially speaking, this condition guarantees that after controlling for
relevant covariates the treatment assignment is completely randomized. Under
this condition, $Q_{d}(\tau;z)$ is identifiable and can be recast as the
solution to the following program:

$\displaystyle
Q_{d}(\tau;\cdot)\in\arg\min_{q(\cdot)}\mathbb{E}\left[\rho_{\tau}(Y-q(X))-\rho_{\tau}(Y)\mid
D=d\right],$ (3)

where $\rho_{\tau}(u)=u(\tau-\mathbf{1}\\{u\leq 0\\})$ is the so-called check-
loss and the minimum is taken over all measurable functions $q(\cdot)$ of $X$
(Koenker,, 2005; Angrist et al.,, 2006). While unconfoundedness of treatment
assignments is a standard condition in the literature on causal inference, it
cannot be verified from the data alone. Rubin, (2009) argues that
unconfoundedness is more plausible when $X$ is a rich set of covariates. This
motivates us to frame our problem as a high-dimensional statistical problem
with predictors $X\in\mathbb{R}^{p}$ whose dimension $p$ exceeds the sample
size $n$.

The convex optimization program (3) poses already a formidable challenge in
low dimensions and to make it tractable in high dimensions we need to impose
further structural constraints:

######  Condition 2 (Sparse linear quantile regression function).

Let $\mathcal{T}$ be a compact subset of $(0,1)$. The CQF of $Y_{d}\mid X=z$
is given by $Q_{d}(\tau;z)=z^{\prime}\theta_{d}(\tau)$ and
$\sup_{\tau\in\mathcal{T}}\left\|\theta_{d}(\tau)\right\|_{0}\ll p\wedge n$.

In principle, this condition can be relaxed to approximate linearity and
approximate sparsity similar to Belloni et al., 2019b , but we do not pursue
the technical refinements in this direction. Under Conditions 1 and 2, the
program (3) reduces to the linear quantile regression program

$\displaystyle\theta_{d}(\tau)\in\arg\min_{\theta\in\mathbb{R}^{p}}\mathbb{E}\left[\rho_{\tau}(Y-X^{\prime}\theta)-\rho_{\tau}(Y)\mid
D=d\right],$ (4)

and the HQTE curve is identified as

$\displaystyle\alpha(\tau;z)=z^{\prime}\theta_{1}(\tau)-z^{\prime}\theta_{0}(\tau).$
(5)

Despite the linearity condition, the HQTE curve in (5) is flexible and can
capture three different aspects of treatment heterogeneity. First, by keeping
$z\in\mathbb{R}^{p}$ fixed and varying only the quantile levels
$\tau\in\mathcal{T}$ we can investigate treatment effect heterogeneity across
different quantile levels. Second, by keeping $\tau\in\mathcal{T}$ fixed and
varying $z\in\mathbb{R}^{p}$ we can analyze individual treatment effects for
individuals characterized by different covariates $z$. Third, by keeping
$\tau\in\mathcal{T}$ fixed and letting $z\in\mathbb{R}^{p}$ be a sparse
contrast we can identify differential effects of treatments in different sub-
populations characterized by a few pre-treatment covariates (e.g. race,
marriage status, gender, socioeconomic status, etc.).

##  3 Methodology

In this section, we introduce the rank-score debiasing procedure for
estimating the HQTE curve. We show that the estimator solves a bias-variance
trade-off problem and discuss its relation to Neyman orthogonalization
(Neyman,, 1959; Belloni et al., 2019b, ).

###  3.1 The rank-score debiasing procedure

Let $\\{(Y_{i},D_{i},X_{i})\\}_{i=1}^{n}$ be a random sample of response
variable $Y$, treatment indicator $D$, and covariates $X$. Denote by
$f_{Y_{d}|X}$ the conditional density of $Y_{d}\mid X$, $d\in\\{0,1\\}$. To
simplify notation, write
$f_{i}(\tau)=f_{Y_{D_{i}}|X}(X_{i}^{\prime}\theta_{D_{i}}(\tau)|X_{i})$,
$i=1,\ldots,n$. Moreover, assume that the first $n_{0}$ observations belong to
the control group and the remaining $n_{1}=n-n_{0}$ observations to the
treatment group.

Step 1. For $d\in\\{0,1\\}$, compute pilot estimates of $\theta_{d}(\tau)$ as
the solution of the $\ell_{1}$-penalized quantile regression program,

$\displaystyle\hat{\theta}_{d}(\tau)\in\arg\min_{\theta\in\mathbb{R}^{p}}\left\\{\sum_{i:D_{i}=d}\rho_{\tau}(Y_{i}-X_{i}^{\prime}\theta)+\lambda_{d}\|\theta\|_{1}\right\\},$
(6)

where $\lambda_{d}>0$ is a regularization parameter. Use the pilot estimates
$\hat{\theta}_{d}(\tau)$ to estimate the conditional densities $f_{i}(\tau)$
as

$\displaystyle\hat{f}_{i}(\tau):=\begin{cases}\frac{2h}{X_{i}^{\prime}\hat{\theta}_{1}(\tau+h)-X_{i}^{\prime}\hat{\theta}_{1}(\tau-h)},&i\in\\{j:D_{j}=1\\}\\\
\frac{2h}{X_{i}^{\prime}\hat{\theta}_{0}(\tau+h)-X_{i}^{\prime}\hat{\theta}_{0}(\tau-h)},&i\in\\{j:D_{j}=0\\},\end{cases}$
(7)

where $h>0$ is a bandwidth parameter. We discuss the choice of $\lambda_{d}$
and $h$ in Sections 5.1 and 5.3.

Step 2. Solve the _rank-score debiasing program_ with plug-in estimates of the
conditional densities from Step 1,

$\displaystyle\widehat{w}(\tau;z)\in\arg\min_{w\in\mathbb{R}^{n}}\left\\{\sum_{i=1}^{n}w_{i}^{2}\hat{f}_{i}^{-2}(\tau):\>\left\|z-\frac{1}{\sqrt{n}}\sum_{i:D_{i}=d}w_{i}X_{i}\right\|_{\infty}\leq\frac{\gamma_{d}}{n},\>\>d\in\\{0,1\\}\right\\},$
(8)

where the $\gamma_{d}>0$ are tuning parameters. We discuss the choice of
$\gamma_{d}$ in Section 5.2.

Step 3. Define the _rank-score debiased estimator of the CQF_ as

$\displaystyle\widehat{Q}_{d}(\tau;z):=z^{\prime}\hat{\theta}_{d}(\tau)+\frac{1}{\sqrt{n}}\sum_{i:D_{i}=d}\widehat{w}_{i}(\tau;z)\hat{f}_{i}^{-1}(\tau)\big{(}\tau-\mathbf{1}\\{Y_{i}\leq
X_{i}^{\prime}\hat{\theta}_{d}(\tau)\\}\big{)}.$ (9)

Step 4. Define the _rank-score debiased estimator of the HQTE curve_ as

$\displaystyle\widehat{\alpha}(\tau;z):=\widehat{Q}_{1}(\tau;z)-\widehat{Q}_{0}(\tau;z).$

and construct an asymptotic $95\%$ confidence interval of $\alpha(\tau;z)$ as

$\displaystyle\left[\widehat{\alpha}(\tau;z)\pm
1.96\times\sqrt{\frac{\tau(1-\tau)}{n}\sum_{i=1}^{n}\widehat{w}_{i}^{2}(\tau;z)\hat{f}_{i}^{-2}(\tau)}\>\right].$

Steps 2 and 3 constitute the core of the rank-score debiasing procedure. In
Step 2 we compute quantile-specific debiasing weights and in Step 3 we augment
the estimated conditional quantile function $z^{\prime}\hat{\theta}_{d}(\tau)$
with a bias correction based on these weights. This bias correction addresses
the penalization bias in $z^{\prime}\hat{\theta}_{d}(\tau)$, because the
$\ell_{1}$-penalty introduces a regularization bias by shrinking coefficients
in $\hat{\theta}_{d}(\tau)$ towards zero. Also, since the quantile regression
vector $\hat{\theta}_{d}(\tau)$ is based on the observed covariates
$\\{X_{i}:D_{i}=d\\}$ alone, estimating $Q_{d}(\tau;z)$ as
$z^{\prime}\hat{\theta}_{d}(\tau)$ introduces a sort of mismatch bias. The
more $z$ differs from a typical covariate in $\\{X_{i}:D_{i}=d\\}$ the larger
is this bias. We refer to our estimator as the rank-score debiased estimator,
because its key component is a weighted sum of quantile regression rank scores
with weights that approximately match the covariates.

###  3.2 Heuristic explanation in terms of a bias-variance trade-off

The rank-score debiased estimator can be motivated in terms of a bias-variance
trade-off. This perspective offers a first glimpse at its theoretical
properties.

Let $\theta\in\mathbb{R}^{p}$ and $w\in\mathbb{R}^{n}$ be arbitrary. To
simplify notation, write
$f_{i}(\tau)=f_{Y_{D_{i}}|X}(X_{i}^{\prime}\theta_{D_{i}}(\tau)|X_{i})$ and
$F_{i}(\tau)=F_{Y_{D_{i}}|X}(X_{i}^{\prime}\theta|X_{i})$ for $i=1,\ldots,n$.
Define $\varphi_{i}(\theta)=\mathbf{1}\\{Y_{i}\leq
X_{i}^{\prime}\theta\\}-\mathbf{1}\\{Y_{i}\leq
X_{i}^{\prime}\theta_{d}(\tau)\\}$ and note that
$\mathbb{E}[f_{i}^{-1}(\tau)\varphi_{i}(\theta)\mid
X_{i}]=f_{i}^{-1}(\tau)\big{(}F_{Y_{D_{i}}|X}(X_{i}^{\prime}\theta|X_{i})-F_{i}(\tau)\big{)}$.
Thus, a first-order Taylor approximation at $\theta=\theta_{D_{i}}(\tau)$
yields

$\displaystyle\frac{1}{\sqrt{n}}\sum_{i:D_{i}=d}w_{i}\mathbb{E}\left[f_{i}^{-1}(\tau)\varphi_{i}(\tau)\mid
X_{i}\right]=\frac{1}{\sqrt{n}}\sum_{i:D_{i}=d}w_{i}X_{i}^{\prime}\big{(}\theta-\theta_{d}(\tau)\big{)}+a_{n}(\theta),$

where
$|a_{n}(\theta)|\leq\big{\|}\frac{1}{\sqrt{n}}\sum_{i:D_{i}=d}w_{i}f_{i}^{-1}(\tau)\xi_{i,\tau}X_{i}X_{i}\big{\|}_{op}\big{\|}\theta-\theta_{d}(\tau)\big{\|}_{2}^{2}$
and $\xi_{i,\tau}=f_{Y_{D_{i}}|X}^{\prime}(X_{i}^{\prime}\xi|X_{i})$ with
$\xi$ an intermediate point between $\theta$ and $\theta_{D_{i}}(\tau)$.
Suppose that this identity remains (approximately) true for
$\theta=\hat{\theta}_{d}(\tau)$. Then, re-arranging this expansion leads to

$\displaystyle\begin{split}&z^{\prime}\hat{\theta}_{d}(\tau)+\frac{1}{\sqrt{n}}\sum_{i:D_{i}=d}w_{i}f_{i}^{-1}(\tau)\big{(}\tau-\mathbf{1}\\{Y_{i}\leq
X_{i}^{\prime}\hat{\theta}_{d}(\tau)\\}\big{)}\\\
&\quad{}\quad{}\quad{}=z^{\prime}\theta_{d}(\tau)+\frac{1}{\sqrt{n}}\sum_{i:D_{i}=d}w_{i}f_{i}^{-1}(\tau)\big{(}\tau-\mathbf{1}\\{Y_{i}\leq
X_{i}^{\prime}\theta_{d}(\tau)\\}\big{)}\\\
&\quad{}\quad{}\quad{}\quad{}\quad{}\quad{}+\left(z-\frac{1}{\sqrt{n}}\sum_{i:D_{i}=d}w_{i}X_{i}\right)^{\prime}\big{(}\hat{\theta}_{d}(\tau)-\theta_{d}(\tau)\big{)}+a_{n}\big{(}\hat{\theta}_{d}(\tau)\big{)}+b_{n}\big{(}\hat{\theta}_{d}(\tau)\big{)},\end{split}$
(10)

where
$b_{n}(\theta)=-\frac{1}{\sqrt{n}}\sum_{i:D_{i}=d}w_{i}\big{(}f_{i}^{-1}(\tau)\varphi_{i}(\theta)-\mathbb{E}[f_{i}^{-1}(\tau)\varphi_{i}(\theta)\mid
X_{i}]\big{)}$.

If $\hat{\theta}_{d}(\tau)$ is consistent for $\theta_{d}(\tau)$ and if the
remainder terms $a_{n}\big{(}\hat{\theta}_{d}(\tau)\big{)}$ and
$b_{n}\big{(}\hat{\theta}_{d}(\tau)\big{)}$ can be shown to be asymptotically
negligible, then the statistical behavior of the left hand side of eq. (10) is
governed by the first three terms on the right hand side. In particular, the
first term on the right hand side, $z^{\prime}\theta_{d}(\tau)$, is
deterministic, the second term has mean zero and variance
$\tau(1-\tau)n^{-1}\sum_{i:D_{i}=d}w_{i}^{2}f_{i}^{-2}(\tau)$ (expectations
taken conditionally on the $X_{i}$’s), and the third term can be upper bounded
by
$\big{\|}z-\frac{1}{\sqrt{n}}\sum_{i:D_{i}=d}w_{i}X_{i}\big{\|}_{\infty}\big{\|}\hat{\theta}_{d}(\tau)-\theta_{d}(\tau)\big{\|}_{1}$.
Since the weights $w$ are arbitrary, we can choose them to fine-tune the
statistical behavior of the left hand side of eq. (10). Given above
observations, it is natural to seek weights $w$ that minimize the variance
$\tau(1-\tau)n^{-1}\sum_{i:D_{i}=d}w_{i}^{2}f_{i}^{-2}(\tau)$ while
controlling the bias term
$\big{\|}z-\frac{1}{\sqrt{n}}\sum_{i:D_{i}=d}w_{i}X_{i}\big{\|}_{\infty}$. The
rank-score debiasing program (8) with plug-in estimates $\hat{f}_{i}(\tau)$
can be viewed as a feasible sample version of this constrained minimization
problem. Since the weights are chosen to minimize the variance of the right
hand side, we expect that the rank-score balanced estimator can be more
efficient than other debiasing procedures. We emphasize that the theoretical
analysis the rank-score debiased estimator does not rely on this Taylor
expansion because it is impossible to bound the remainder terms uniformly in
$w\in\mathbb{R}^{n}$ as $n$ diverges.

###  3.3 Connection to Neyman orthogonalization

Our algorithm can also be rationalized as an approximate Neyman
orthogonalization procedure (Neyman,, 1959; Belloni et al., 2019b, ;
Chernozhukov et al.,, 2018).

Given the target $Q_{d}(\tau;z)$, one may interpret the true quantile
regression coefficient $\theta_{d}(\tau)$ as a nuisance parameter, say
$\eta_{0}\equiv\theta_{d}(\tau)$. To carry out valid inference on
$Q_{d}(\tau;z)$ when the high-dimensional nuisance parameter $\eta_{0}$ cannot
be estimated at $\sqrt{n}$-rate, one then seeks a score function
$\psi(q,\eta)$ such that for all $\eta$ in a (shrinking) neighborhood
$\mathcal{N}_{n}$ of $\eta_{0}$ and a null sequence $(\delta_{n})_{n\geq 1}$,

$\displaystyle\mathbb{E}\big{[}\psi(Q_{d}(\tau;z),\eta_{0})\mid
X\big{]}=0\quad{}\quad{}\mathrm{and}\quad{}\quad{}\sup_{\eta\in\mathcal{N}_{n}}\left|\frac{\partial}{\partial\eta}\mathbb{E}\big{[}\psi(Q_{d}(\tau;z),\eta_{0})\mid
X\big{]}(\eta-\eta_{0})\right|\leq\delta_{n}n^{-1/2}.$ (11)

These equations are known as Neyman near-orthogonality conditions
(Chernozhukov et al.,, 2018, Section 3.2). The neighborhood $\mathcal{N}_{n}$
is also called the nuisance realization set and chosen such that it contains
the estimated nuisance parameter $\hat{\eta}$ with high probability. Given the
definition of the rank-score debiased estimator in eq. (9), a natural choice
for the score function is

$\displaystyle\psi_{w}(q,\eta):=q-z^{\prime}\eta-\frac{1}{\sqrt{n}}\sum_{i:D_{i}=d}w_{i}f_{i}^{-1}(\tau)(\tau-\mathbf{1}\big{(}\\{Y_{i}\leq
X_{i}^{\prime}\eta\\}\big{)},$

where $w\in\mathbb{R}^{n}$ is a tuning parameter to be chosen later. One
easily verifies that the score function $\psi_{w}$ satisfies the first
equality in eq. (11) for all $w\in\mathbb{R}^{n}$. Furthermore, provided that
$w\in\mathbb{R}^{n}$ satisfies the box-constraint in program (8) and that the
nuisance realization set can be chosen as
$\mathcal{N}_{n}=\big{\\{}\eta\in\mathbb{R}^{p}:\|\eta-\eta_{0}\|_{1}\leq\frac{\delta_{n}}{\gamma_{d}}n^{1/2}\big{\\}}$,
the second inequality in eq. (11) holds as well: Indeed, for all
$\eta\in\mathcal{N}_{n}$, by Hölder’s inequality,

$\displaystyle\left|\frac{\partial}{\partial\eta}\mathbb{E}\big{[}\psi_{w}(Q_{d}(\tau;z),\eta_{0})\mid
X\big{]}(\eta-\eta_{0})\right|=\left|\left(z-\frac{1}{\sqrt{n}}\sum_{i:D_{i}=d}w_{i}X_{i}\right)^{\prime}(\eta-\eta_{0})\right|\leq\delta_{n}n^{-1/2}.$

Next, denote by $\widehat{Q}_{d}(\tau;z,w)$ the generalized method of moment
estimator that solves $\sum_{i:D_{i}=d}\psi_{w}\big{(}\\\
\widehat{Q}_{d}(\tau;z,w),\hat{\eta}\big{)}=0$. Conditionally on the
$X_{i}$’s, $\widehat{Q}_{d}(\tau;z,w)$ has asymptotic variance
$\tau(1-\tau)n^{-1}\sum_{i:D_{i}=d}\\\ w_{i}^{2}f_{i}^{-2}(\tau)$ (e.g.
Chernozhukov et al.,, 2018, Section 3.2). Since $w\in\mathbb{R}^{n}$ is
arbitrary, it is sensible to choose $w$ to minimize this asymptotic variance.
Hence, the rank-score debiasing algorithm with plug-in estimates
$\hat{f}_{i}(\tau)$ can be viewed as a feasible sample version of this
approximate Neyman orthogonalization procedure. Intuitively, the box-
constraint in program (8) relaxes the strict Neyman orthogonality condition
since in high dimensions one can not hope to match $z$ exactly with a linear
combination of the $X_{i}$’s. Furthermore, the inverse-density weighting of
the weights in the expression $\sum_{i:D_{i}=d}w_{i}^{2}f_{i}^{-2}(\tau)$
ensures that observations associated with low density at the $\tau$th quantile
are given smaller debiasing weights.

##  4 Theoretical analysis

In this section we establish joint asymptotic normality of the HQTE process,
propose consistent estimators of its asymptotic covariance function, and
discuss the duality theory of the rank-score debiasing program which underlies
the theoretical results.

###  4.1 Regularity conditions

Throughout, we assume that $\\{(Y_{i},D_{i},X_{i})\\}_{i=1}^{n}$ are i.i.d.
copies of $(Y,D,X)$. Recall that $Y=DY_{1}+(1-D)Y_{0}\in\mathbb{R}$, where
$Y_{1}$ and $Y_{0}$ are potential outcomes, $D\in\\{0,1\\}$, and
$X\in\mathbb{R}^{p}$. For examples of quantile regression models that satisfy
below conditions, we refer to Section 4.2.

######  Condition 3 (Sub-Gaussian predictors).

$X\in\mathbb{R}^{p}$ is a sub-Gaussian vector, i.e.
$\|X-\mathbb{E}[X]\|_{\psi_{2}}\lesssim\big{(}\mathbb{E}[(X^{\prime}u)^{2}]\big{)}^{1/2}$
for all $u\in\mathbb{R}$.

Condition 3 is standard in high-dimensional statistics. We introduce it to
analyze the rank-score debiasing program (8), but it also simplifies the
theoretical analysis of the quantile regression program (6). The specific
formulation of sub-Gaussianity is convenient because it allows us to relate
higher moments of (sparse) linear combinations $X^{\prime}u$ to (sparse)
eigenvalues of their covariance and second moment matrix (i.e. design matrix).

We require the following conditions on the conditional quantiles and density
of $Y_{d}$ given $X$:

######  Condition 4 (Sparsity and Lipschitz continuity of
$\tau\mapsto\theta_{d}(\tau)$).

Let $\mathcal{T}$ be compact subset of $(0,1)$.

  * (i)

There exists $s_{\theta}\geq 1$ such that
$\sup_{d\in\\{0,1\\}}\sup_{\tau\in\mathcal{T}}\big{|}T_{\theta_{d}}(\tau)\big{|}\leq
s_{\theta}$ for
$T_{\theta_{d}}(\tau)=\mathrm{support}\big{(}\theta_{d}(\tau)\big{)}$;

  * (ii)

There exists $L_{\theta}\geq 1$ such that
$\sup_{d\in\\{0,1\\}}\|\theta_{d}(\tau)-\theta_{d}(\tau^{\prime})\|_{2}\leq
L_{\theta}|\tau-\tau^{\prime}|$ for all $\tau,\tau^{\prime}\in\mathcal{T}$.

######  Condition 5 (Boundedness and Lipschitz continuity of $f_{Y_{d}|X}$).

Let $a,b,x\in\mathbb{R}^{p}$ be arbitrary.

  * (i)

There exists $\bar{f}\geq 1$ such that
$\sup_{d\in\\{0,1\\}}f_{Y_{d}|X}(a|x)\leq\bar{f}$;

  * (ii)

There exists $\underline{f}>0$ such that
$\inf_{d\in\\{0,1\\}}\inf_{\tau\in\mathcal{T}}f_{Y_{d}|X}(x^{\prime}\theta_{d}(\tau)|x)\geq\underline{f}$;

  * (iii)

There exists $L_{f}\geq 1$ such that
$\sup_{d\in\\{0,1\\}}\left|f_{Y_{d}|X}(x^{\prime}a|x)-f_{Y_{d}|X}(x^{\prime}b|x)\right|\leq
L_{f}|x^{\prime}a-x^{\prime}b|$.

######  Condition 6 (Differentiability of $\tau\mapsto Q_{d}(\tau;X)$).

Let $\mathcal{T}$ be a compact subset of $(0,1)$. The CQF $Q_{d}(\tau;X)$ is
three times boundedly differentiable on $\mathcal{T}$, i.e. there exists
$C_{Q}\geq 1$ such that
$\sup_{d\in\\{0,1\\}}\left|Q_{d}^{\prime\prime\prime}(\tau;x)\right|\leq
C_{Q}$ for all $x\in\mathbb{R}^{p}$ and $\tau\in\mathcal{T}$.

Conditions 4 and 5 are common in the literature on high-dimensional quantile
regression (Belloni and Chernozhukov,, 2011; Chao et al.,, 2017; Belloni et
al., 2019b, ; Wang and He,, 2021). They are relevant for establishing weak
convergence of the rank-score debiased HQTE process to a Gaussian process in
$\ell^{\infty}(\mathcal{T})$. Conditions 5 $(i)$ and $(ii)$ are only needed
for the theoretical analysis of program (8) and for establishing uniform (in
$\tau\in\mathcal{T}$) consistency of the non-parametric estimates of the
conditional densities in (7); for all other purposes they can be dropped. If
one is only interested in consistency and asymptotic normality of a single (or
finitely many) quantile level(s), one can also drop Conditions 4 $(ii)$ and
Condition 5 $(iii)$. Condition 6 was introduced recently in Belloni et al.,
2019b  as part of the sufficient conditions for establishing consistency of
the non-parametric estimates of the conditional densities in (7). It might be
possible to relax this condition to $Q_{d}(\tau;x)$ belonging to a Hölder
class of functions, which is a common assumption in non-parametric (quantile)
spline estimation (He and Shi,, 1994; He et al.,, 2013).

The next two definitions and conditions are variations of canonical
assumptions for high-dimensional regression models.

######  Definition 1 ($s$-sparse maximum eigenvalues).

We define the $s$-sparse maximum eigenvalues of the population and sample
design matrices by

$\displaystyle\varphi_{\max,d}(s):=\sup_{u:\|u\|_{0}\leq
s}\frac{\mathbb{E}[(X^{\prime}u)^{2}\mathbf{1}\\{D=d\\}]}{\|u\|_{2}^{2}}\hskip
20.0pt\mathrm{and}\hskip
20.0pt\widehat{\varphi}_{\max,d}(s):=\sup_{u:\|u\|_{0}\leq
s}\frac{n^{-1}\sum_{i:D_{i}=d}(X_{i}^{\prime}u)^{2}}{\|u\|_{2}^{2}}.$

######  Condition 7 (Bounds on maximum eigenvalues).

There exists an absolute constant $\varphi_{\max}\geq 1$ such that

$\displaystyle\varphi_{\max,d}\big{(}n_{d}/\log(n_{d}p)\big{)}\vee\widehat{\varphi}_{\max,d}\big{(}n_{d}/\log(n_{d}p)\big{)}\leq\varphi_{\mathrm{max}},\quad{}\quad{}d\in\\{0,1\\}.$

Under Condition 3 and for $\log p=o(n_{d})$ one can upper bound the empirical
maximal eigenvalue
$\widehat{\varphi}_{\max,d}\big{(}n_{d}/\log(n_{d}p)\big{)}$ by a constant
multiple of $\varphi_{\max,d}\big{(}n_{d}/\log(n_{d}p)\big{)}$ with
probability tending to 1 (e.g. apply Lemma 14 in the Appendix). Hence,
Condition 7 is first and foremost a condition on the maximum eigenvalue of the
population design matrix.

To state the next definition recall that for $J\subseteq\\{1,\ldots,p\\}$,
$q\geq 1$, and $\vartheta\in[0,\infty]$ the cone of $(J,\vartheta)$-dominant
coordinate is defined as
$C^{p}_{q}(J,\vartheta):=\left\\{u\in\mathbb{R}^{p}:\|u_{J^{c}}\|_{q}\leq\vartheta\|u_{J}\|_{q}\right\\}$.

######  Definition 2 ($(\omega,\vartheta,\varrho)$-restricted minimum
eigenvalue of the design matrix).

Let $\mathcal{T}$ be a compact subset of $(0,1)$ and
$\omega,\vartheta,\varrho\geq 0$. We define the
$(\omega,\vartheta,\varrho)$-restricted minimum eigenvalue of the design
matrix as

$\displaystyle\kappa_{\omega}(\vartheta,\varrho):=\min_{d\in\\{0,1\\}}\inf_{\tau\in\mathcal{T}}\inf_{\|\zeta\|_{2}\leq\varrho\>}\inf_{u\in
C^{p}_{1}(T_{\theta}(\tau),\vartheta)\cap\partial
B^{p}_{2}(0,1)}\mathbb{E}\left[f_{Y|X}^{\omega}\big{(}X^{\prime}\theta_{0}(\tau)+X^{\prime}\zeta|X\big{)}(X^{\prime}u)^{2}\mathbf{1}\\{D=d\\}\right].$

To simplify notation we write
$\kappa_{\omega}(\vartheta):=\kappa_{\omega}(\vartheta,0)$.

######  Condition 8 ($\varrho_{n}$-restricted identifiability of
$\theta_{d}(\tau)$).

Let $\mathcal{T}$ be a compact subset of $(0,1)$ and $(\varrho_{n})_{n\geq 1}$
a null sequence. The quantile regression vectors $\theta_{d}(\tau)$ with
$d\in\\{0,1\\}$ and $\tau\in\mathcal{T}$ are $\varrho_{n}$-restricted
identifiable, i.e. $\kappa_{1}(2)>0\ \mathrm{and}\
\kappa_{1}(2,\varrho_{n})\gtrsim\kappa_{1}(2).$

Condition 8 guarantees that the objective function of the $\ell_{1}$-penalized
quantile regression program (6) can be locally minorized by a quadratic
function. To the best of our knowledge this identifiability condition for
high-dimensional quantile regression vectors is new. We use it with
$\varrho_{n}\asymp\sqrt{s_{\theta}(\log np)/n}$. For this choice of
$\varrho_{n}$, Condition 8 is milder than the restricted identifiability and
nonlinearity condition D.5 in Belloni and Chernozhukov, (2011) and also
slightly less restrictive than Condition (C1) in Wang and He, (2021). For a
comparison of these conditions, see Remark 1 in Wang and He, (2021) and
Section F.2 in the Supplementary Materials.

The last set of definitions and conditions concern the dual of the rank-score
debiasing program 8. Readers may skip over these conditions and return to them
after having read Section 4.4.

######  Definition 3 ($\epsilon$-approximation).

Let $\epsilon\geq 0$. We call a vector $\tilde{v}\in\mathbb{R}^{p}$ an
$\epsilon$-approximation of $v\in\mathbb{R}^{p}$ if
$\|v-\tilde{v}\|_{2}\leq\epsilon\|v\|_{2}$.

######  Condition 9 (Sparse $\epsilon_{n}$-approximate solution to the
population dual).

For $z\in\mathbb{R}^{p}$, $\tau\in\mathcal{T}$ and $d\in\\{0,1\\}$ define

$\displaystyle
v_{d}(\tau;z):=-2\mathbb{E}\left[f_{Y_{d}|X}^{2}\big{(}X^{\prime}\theta_{d}(\tau)|X\big{)}XX^{\prime}\mathbf{1}\\{D=d\\}\right]^{-1}z.$

Let $(\epsilon_{n})_{n\geq 1}$ be a null sequence and
$\big{(}\tilde{v}_{d,n}(\tau;z)\big{)}_{n\geq 1}$ the associated collection of
$\epsilon_{n}$-approximations of $v_{d}(\tau;z)$. We assume that there exists
$(s_{v,n})_{n\geq 1}$ such that

$\displaystyle\sup_{d\in\\{0,1\\}}\sup_{\tau\in\mathcal{T}}\big{|}T_{v_{d,n}}(\tau)\big{|}\leq
s_{v,n}\ll n\wedge
p,\quad{}\quad{}\mathrm{where}\quad{}\quad{}T_{v_{d,n}}(\tau)=\mathrm{support}\big{(}\tilde{v}_{d,n}(\tau;z)\big{)}.$

We drop the subscript $n$ on $s_{v,n}$ and $v_{d,n}(\tau;z)$ if this does not
cause confusion.

Condition 9 is a technical condition that allows us to analyze the rank-score
debiasing weights. The plausibility of Condition 9 depends crucially on the
choice of $(\epsilon_{n})_{n\geq 1}$. Intuitively, the larger
$\epsilon_{n}\geq 0$, the easier it is to find a sparse
$\epsilon_{n}$-approximation $\tilde{v}_{d}(\tau;z)$ of $v_{d}(\tau;z)$.
Indeed, if $\epsilon_{n}\geq 1$, then one may take
$\tilde{v}_{d}(\tau;z)\equiv 0$ with $s_{v}=0$. In contrast, if
$\epsilon_{n}=0$, then, necessarily, $\tilde{v}_{d}(\tau;z)=v_{d}(\tau;z)$ and
$s_{v}=\|v_{d}(\tau;z)\|_{0}$, which may or may not be less than $n\wedge p$.
Our theoretical results hold for any null sequence $\epsilon_{n}\lesssim
1/\sqrt{s_{v}}$. Typically, we choose $s_{v}\asymp\log n$, and, hence,
Definition 3 and Condition 9 combine the notion of sieve estimators from
classical statistics (e.g. Chen,, 2007) with the concept of compressibility
from the literature on compressive sensing (e.g. Foucart and Rauhut,, 2013).
We provide concrete examples and high-level conditions under which Condition 9
holds in Section 4.2. To simplify the presentation, above definition and
condition are stated somewhat informal. The rigorous formulations can be found
in Section G.2 in the Supplementary Materials.

######  Condition 10 (Identifiability of $\tilde{v}_{d}(\tau;z)$).

The sparse $\epsilon_{n}$-approximate solution to the population dual
$\tilde{v}_{d}(\tau;z)$ is identifiable, i.e. $\kappa_{2}(\infty)>0$.

Condition 10 guarantees that the objective function of the dual of program (8)
can be locally minorized by a quadratic function.

###  4.2 Examples of simple sufficient conditions

We illustrate the general Conditions 3–10 with some simple sufficient
conditions. We emphasize that the conditions of Section 4.1 are significantly
more general than the examples discussed here.

######  Example 1 (Location model with Gaussian predictors and autoregressive
covariance structure).

Consider the location model

$\displaystyle
Y_{d}=\alpha_{d}+X^{\prime}\beta_{d}+\varepsilon,\quad{}\quad{}X\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\varepsilon,\quad{}\quad{}d\in\\{0,1\\},$

where $\varepsilon\sim N(0,\sigma_{\varepsilon}^{2})$,
$\sigma_{\varepsilon}>0$ fixed, $X\sim N(0,\Sigma)$, and smallest and largest
eigenvalues of $\Sigma\in\mathbb{R}^{p\times p}$ bounded from below by
$\underline{\kappa}>0$ and from above by $\bar{\varphi}<\infty$. Moreover,
suppose that the precision matrix
$\Sigma^{-1}\equiv\Omega=(\omega_{jk})_{j,k=1}^{p}$ has bandwidth $1\leq q<p$,
i.e. $\omega_{jk}=0$ if $k<j-q$ or $k>j+q$.

######  Lemma 1.

Let $z\in\mathbb{R}^{p}$ be sparse with $\|z\|_{0}\leq s_{z}$ and
$\mathcal{T}=[\xi,1-\xi]$. Under the design in Example 1, Condition 3–10 are
satisfied with

$\displaystyle
s_{\theta}\leq\max_{d\in\\{0,1\\}}\|\beta_{d}\|_{0}+1,\quad{}L_{\theta}=\sigma_{\varepsilon}/\xi\vee
1,\quad{}\bar{f}=1/\sqrt{2\pi\sigma^{2}_{\varepsilon}}\vee
1,\quad{}\underline{f}=\sqrt{\xi}/\sqrt{2\pi\sigma^{2}_{\varepsilon}},\quad{}s_{v}\leq(q+1)s_{z},$
$\displaystyle L_{f}=\sqrt{e/(2\pi\sigma^{4}_{\varepsilon})}\vee
1,\quad{}C_{Q}=4\sigma_{\varepsilon}/\xi^{4},\quad{}\varphi_{\max}=\bar{\varphi},\quad{}\kappa_{1}(2)\geq\sqrt{\xi}\underline{\kappa}/\sqrt{2\pi\sigma^{2}_{\varepsilon}},\quad{}\kappa_{2}(\infty)\geq\xi\underline{\kappa}/(2\pi\sigma^{2}_{\varepsilon}),$
$\displaystyle\varrho_{n}=o(1),\quad{}\epsilon_{n}=0.$

In Example 1 the covariance structure and the sparsity of $z$ guarantee that
$v_{d}(\tau;z)$ is sparse. Hence, Condition 9 is trivially satisfied. In the
next two examples we only require $v_{d}(\tau;z)$ to lie in some cone of
dominant coordinates. This is a mild assumption and allows $v_{d}(\tau;z)$ to
be dense and/ or weakly sparse (see also Lemma 4 below).

######  Example 2 (Location model with Gaussian predictors).

Consider the location model

$\displaystyle
Y_{d}=\alpha_{d}+X^{\prime}\beta_{d}+\varepsilon,\quad{}\quad{}X\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\varepsilon,\quad{}\quad{}d\in\\{0,1\\},$

where $\varepsilon\sim N(0,\sigma_{\varepsilon}^{2})$,
$\sigma_{\varepsilon}>0$ fixed, $X\sim N(0,\Sigma)$, and smallest and largest
eigenvalues of $\Sigma\in\mathbb{R}^{p\times p}$ bounded from below by
$\underline{\kappa}>0$ and from above by $\bar{\varphi}<\infty$.

######  Lemma 2.

Let $\mathcal{T}=[\xi,1-\xi]$, $c_{0}\in(0,\infty]$, and
$J\subseteq\\{1,\ldots,p\\}$ with $|J|\leq s$. Suppose that $v_{d}(\tau;z)\in
C^{p}_{1}(J,c_{0})$ for $d\in\\{0,1\\}$ and $\tau\in\mathcal{T}$. Under the
design in Example 2, Condition 3–10 are satisfied with

$\displaystyle
s_{\theta}\leq\max_{d\in\\{0,1\\}}\|\beta_{d}\|_{0}+1,\quad{}L_{\theta}=\sigma_{\varepsilon}/\xi\vee
1,\quad{}\bar{f}=1/\sqrt{2\pi\sigma^{2}_{\varepsilon}}\vee
1,\quad{}\underline{f}=\sqrt{\xi}/\sqrt{2\pi\sigma^{2}_{\varepsilon}},\quad{}s_{v}=s\log
n$ $\displaystyle L_{f}=\sqrt{e/(2\pi\sigma^{4}_{\varepsilon})}\vee
1,\quad{}C_{Q}=4\sigma_{\varepsilon}/\xi^{4},\quad{}\varphi_{\max}=\bar{\varphi},\quad{}\kappa_{1}(2)\geq\sqrt{\xi}\underline{\kappa}/\sqrt{2\pi\sigma^{2}_{\varepsilon}},\quad{}\kappa_{2}(\infty)\geq\xi\underline{\kappa}/(2\pi\sigma^{2}_{\varepsilon}),$
$\displaystyle\varrho_{n}=o(1),\quad{}\epsilon_{n}=o\left(c_{0}/\sqrt{\log
n}\right).$

######  Example 3 (Location-scale model with bounded predictors).

Consider the location-scale model

$\displaystyle Y_{d}=X^{\prime}\beta_{d}+\varepsilon\cdot
X^{\prime}\eta_{d},\quad{}\quad{}X\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\varepsilon,\quad{}\quad{}d\in\\{0,1\\},$

where $\varepsilon\sim F$ with twice boundedly differentiable density $f$.
Suppose that the smallest and largest eigenvalues of
$\mathbb{E}[XX^{\prime}]\in\mathbb{R}^{p\times p}$ are bounded from below by
$\underline{\kappa}>0$ and from above by $\bar{\varphi}<\infty$. Furthermore,
suppose that there exit absolute constants $K,\upsilon,\Upsilon>0$ such that
$\max_{1\leq k\leq p}|x^{(k)}|\leq K$ and $0<\upsilon\leq
x^{\prime}\eta\leq\Upsilon<\infty$ for all
$x=(x^{(1)},\ldots,x^{(p)})^{\prime}$ in the range of $X$.

######  Lemma 3.

Let $\mathcal{T}=[\xi,1-\xi]$, $c_{0}\in(0,\infty]$, and
$J\subseteq\\{1,\ldots,p\\}$ with $|J|\leq s$. Suppose that $v_{d}(\tau;z)\in
C^{p}_{1}(J,c_{0})$ for $d\in\\{0,1\\}$ and $\tau\in\mathcal{T}$. Under the
design in Example 3, Condition 3–10 are satisfied with

$\displaystyle
s_{\theta}\leq\max_{d}\|\beta_{d}\|_{0}+\|\eta_{d}\|_{0},\quad{}L_{\theta}=\max_{d}\|\eta_{d}\|_{2}\underline{f},\quad{}\bar{f}=\max_{y}f(y)/\upsilon\vee
1,$
$\displaystyle\underline{f}=\min_{\tau\in\mathcal{T}}f(F^{-1}(\tau))/\Upsilon,\quad{}s_{v}=s\log
n,\quad{}L_{f}=\max_{f}f^{\prime}(y)/\upsilon^{2}\vee 1,$ $\displaystyle
C_{Q}=\max_{y}\left(f^{\prime\prime}(y)/\underline{f}^{4}+3L_{f}^{2}\upsilon^{4}/\underline{f}^{5}\right)\Upsilon,\quad{}\varphi_{\max}=\bar{\varphi},\quad{}\kappa_{1}(2)\geq\underline{f}\underline{\kappa},$
$\displaystyle\kappa_{2}(\infty)\geq\underline{f}^{2}\underline{\kappa},\quad{}\varrho_{n}=o(1),\quad{}\epsilon_{n}=o\left(c_{0}/\sqrt{\log
n}\right).$

In above three examples we have imposed high-level assumptions on
$v_{0}(\tau;d)$ which guarantee that Condition 9 holds. The next lemma
provides more specific and (to some extent) testable sufficient conditions
under which Condition 9 is met.

######  Lemma 4 (Sufficient conditions for sparse $\epsilon_{n}$-approximate
solutions to the population dual).

To simplify notation, write
$A=[A_{1},\ldots,A_{p}]:=\mathbb{E}\left[f_{Y_{d}|X}^{2}\big{(}X^{\prime}\theta_{d}(\tau)|X\big{)}XX^{\prime}\mathbf{1}\\{D=d\\}\right]^{-1}\in\mathbb{R}^{p\times
p}$. For subsets $S,T\subseteq\\{1,\ldots,p\\}$ let
$A_{S,T}\in\mathbb{R}^{|S|\times|T|}$ be the sub-matrix obtained from $A$ by
deleting all rows in $S^{c}$ and columns in $T^{c}$. Denote by
$\sigma_{\min}(A_{S,T})$ the smallest positive singular value of $A_{S,T}$ and
set $\kappa_{\min}(S,c_{0}):=\inf_{u\in C^{p}_{1}(S,c_{0})}\|A_{S}u\|_{2}$ for
$c_{0}\geq 0$.

  * (i)

If each column of $A$ has at most $q\geq 1$ non-zero entries and
$z\in\mathbb{R}^{p}$ has at most $s_{z}\geq 1$ non-zero entries, then
Condition 9 holds with $s_{v}=qs_{z}$ and $\epsilon_{n}\equiv 0$ for all
$n\geq 1$.

  * (ii)

If each column of $A$ lies in $C^{p}_{1}(J,\vartheta)$ for some
$J\subseteq\\{1,\ldots,p\\}$ and $\vartheta\in[0,\infty)$ and
$z\in\mathbb{R}^{p}$ has at most $1\leq s_{z}<|J|$ non-zero entries, then
Condition 9 holds with $s_{v}=|J|\log n$ and $\epsilon_{n}=O\big{(}\vartheta
K(J,z)/\sqrt{\log n}\big{)}$, where $K(J,z)=\max_{k\in
T_{z}}\sqrt{s_{z}}\left\|A_{J,k}\right\|_{1}/\sigma_{\min}(A_{J,T_{z}})$ and
$T_{z}=\mathrm{support}(z)$.

  * (iii)

If each column of $A$ lies in $C^{p}_{1}(J,\vartheta)$ for some
$J\subseteq\\{1,\ldots,p\\}$ and $\vartheta\in[0,\infty)$ and
$z\in\mathbb{R}^{p}$ satisfies $z_{J}\neq 0$, then Condition 9 holds with
$s_{v}=|J|\log n$ and $\epsilon_{n}=O\big{(}\vartheta K(J,z)/\sqrt{\log
n}\big{)}$, where $K(J,z)=(1+c_{0})\max_{1\leq k\leq
p}\sqrt{|J|}\left\|A_{J,k}\right\|_{1}/\kappa_{\min}(J,c_{0})$ and
$c_{0}=\|z_{J^{c}}\|_{1}/\|z_{J}\|_{1}$.

From this lemma we infer that Condition 9 holds whenever the columns of
$\mathbb{E}\big{[}f_{Y_{d}|X}^{2}\big{(}X^{\prime}\theta_{d}(\tau)|X\big{)}\\\
XX^{\prime}\mathbf{1}\\{D=d\\}\big{]}^{-1}\in\mathbb{R}^{p\times p}$ are
(weakly) sparse. Typically, this is the case if most predictors are only
weakly correlated. Moreover, sparsity of $z\in\mathbb{R}^{p}$ is not
necessary; in particular, by part $(iii)$ $\|z\|_{1}=O(1)$ is sufficient. We
illustrate these facts in the following example:

######  Example 4 (Homoscedastic quantile regression model).

Suppose that $Y_{d}=\alpha_{d}+X^{\prime}\beta_{d}+\varepsilon$ with
$X,D,\varepsilon$ independent of each other for all $d\in\\{0,1\\}$. Let
$Q_{\varepsilon}(\tau)$ be the $\tau$th quantile of the error $\varepsilon$
and $0<\mathbb{P}\\{D=1\\}=\pi_{1}=1-\pi_{0}=1-\mathbb{P}\\{D=0\\}<1$. Then,

$\displaystyle
v_{d}(\tau;z)=-2\pi_{d}^{-1}f_{\varepsilon}(Q_{\varepsilon}(\tau))^{-2}\mathbb{E}[XX^{\prime}]^{-1}z.$

From this expression we easily read off the following:

  * (i)

If $z\in\mathbb{R}^{p}$ has at most $s_{z}\geq 1$ non-zero entries and at
least $p-q\geq 1$ entries in $X\sim N(0,\Sigma)$ are independent or $X$
follows an AR(q) process, then Lemma 4 (i) applies.

  * (ii)

If $z\in\mathbb{R}^{p}$ has at most $1\leq s_{z}<q$ non-zero entries and $X$
follows an MA(q) process, then there exist a set $J\subseteq\\{1,\ldots,p\\}$
with $|J|=q$ and $\vartheta\in[0,\infty)$ such that Lemma 4 (ii) applies.

  * (iii)

If $z\in\mathbb{R}^{p}$ satisfies $\|z\|_{1}=O(1)$ and at least $p-q\geq 1$
entries in $X\sim N(0,\Sigma)$ are independent or $X$ follows an AR(q) or
MA(q) process, then there exist a set $J\subseteq\\{1,\ldots,p\\}$ with
$|J|=q$ and $\vartheta\in[0,\infty)$ such that Lemma 4 (iii) applies.

###  4.3 Weak convergence results

In this section we establish weak convergence of the rank-score debiased CQF
and the HQTE processes,

$\displaystyle\left\\{\sqrt{n}\big{(}\widehat{Q}_{d}(\tau;z)-Q_{d}(\tau;z)\big{)}:\tau\in\mathcal{T}\right\\}\quad{}\quad{}\mathrm{and}\quad{}\quad{}\left\\{\sqrt{n}\big{(}\widehat{\alpha}(\tau;z)-\alpha(\tau;z)\big{)}:\tau\in\mathcal{T}\right\\}.$

The large sample properties of these processes are needed whenever one would
like to conduct inference on the HQTE curve on more than just one quantile at
a time. For example, statistical comparisons of the HQTE across different
quantiles require uniform confidence bands that hold for all quantiles under
consideration. Similarly, testing hypotheses about subsets of quantiles
requires constructing rejection regions that hold across these quantiles. In
both cases, process methods provide a natural way of addressing these
problems. We provide concrete examples below.

To formulate the theoretical results we introduce the following operator:

$\displaystyle
H_{d}^{(n)}(\tau_{1},\tau_{2};z):=v_{d}^{\prime}(\tau_{1};z)\mathbb{E}[f_{Y_{d}|X}(X^{\prime}\theta_{d}(\tau_{1})|X)f_{Y_{d}|X}(X^{\prime}\theta_{d}(\tau_{2})|X)XX^{\prime}\mathbf{1}\\{D=d\\}]v_{d}(\tau_{2};z),$

where
$v_{d}(\tau;z)=-2\big{(}\mathbb{E}[f_{Y_{d}|X}^{2}\big{(}X^{\prime}\theta_{d}(\tau)|X\big{)}XX^{\prime}\mathbf{1}\\{D=d\\}]\big{)}^{-1}z$.
Since the dimension $n$ may grow with the sample size $n$, we make the
dependence of $H_{d}^{(n)}(\tau_{1},\tau_{2};z)$ on $n$ explicit.

The following theorem establishes joint asymptotic normality of the rank-score
balanced CQF process.

######  Theorem 1 (Weak convergence of the rank-score debiased CQF process).

Let $\mathcal{T}$ be a compact subset of $(0,1)$. Suppose that Conditions 1–10
hold with $\varrho_{n}=\sqrt{(s_{v}+s_{\theta})\log(np)/n}$ and
$\epsilon_{n}^{2}=O(\sqrt{n}h^{-1}\varrho_{n}^{2}+h^{2})$. In addition,
suppose that $(s_{v}+s_{\theta})^{3}\log^{3}(np)\log^{3}(n)=o(nh^{3})$,
$h^{2}s_{v}=o(1)$, and $\|z\|_{2}=O(1)$. If
$\lambda_{d}\asymp\sqrt{\varphi_{\max}}\sqrt{n\log(np)}$ and
$\gamma_{d}\asymp\|z\|_{2}\left(h^{-1}s_{\theta}\log(np)+\sqrt{n}h^{2}\right)\sqrt{n}$,
then

$\displaystyle\sqrt{n}\big{(}\widehat{Q}_{d}(\cdot;z)-Q_{d}(\cdot;z)\big{)}\leadsto\mathbb{G}_{d}(\cdot\>;z)\quad{}\mathrm{in}\quad{}\ell^{\infty}(\mathcal{T}),$

where $\mathbb{G}_{d}(\cdot\>;z)$ is a centered Gaussian process with
covariance function $(\tau_{1},\tau_{2})\mapsto
H_{d}(\tau_{1},\tau_{2};z):=\lim_{n\rightarrow\infty}\frac{\tau_{1}\wedge\tau_{2}-\tau_{1}\tau_{2}}{4}H_{d}^{(n)}(\tau_{1},\tau_{2};z)$
provided this limit exists pointwise for all
$\tau_{1},\tau_{2}\in\mathcal{T}$.

######  Remark 1 (On the existence of the covariance function).

It is easy to verify that the limit $H_{d}(\tau_{1},\tau_{2};z)$ is finite for
all $\tau_{1},\tau_{2}\in\mathcal{T}$ whenever Condition 10 holds and
$\|z\|_{2}=O(1)$. However, this alone does not imply existence of the limit,
since $H_{d}^{(n)}(\tau_{1},\tau_{2};z)$ may oscillate with the sample size
$n$. Hence, we impose pointwise convergence of
$H_{d}^{(n)}(\tau_{1},\tau_{2};z)$ for all $\tau_{1},\tau_{2}\in\mathcal{T}$
as an additional assumption. In the context of abstract weak convergence
results for classes of functions that may change with the sample size $n$ this
assumption is standard (e.g. van der Vaart and Wellner,, 1996, ch. 2.11.3); in
the context of high-dimensional quantile regression this assumption also
appears in Chao et al., (2017). If $H_{d}^{(n)}(\tau_{1},\tau_{2};z)$ does not
converge pointwise for all $\tau_{1},\tau_{2}\in\mathcal{T}$ weak process
convergence fails, but we still have asymptotic normality of the studentized
rank-score debiased CQF: Indeed, for all (fixed) $\tau\in\mathcal{T}$, Lemma 7
implies that
$\widehat{\sigma}_{2}^{-1}(\tau;z)\sqrt{n}\big{(}\widehat{Q}_{d}(\tau;z)-Q_{d}(\tau;z)\big{)}\leadsto
N(0,1)$, where $\widehat{\sigma}_{2}(\tau;z)$ is defined in eq. (16).

Assume, for a moment, that dimension $p$ is fixed. Then, Theorem 1 implies
that

$\displaystyle\sqrt{n}\big{(}\widehat{Q}_{d}(\tau;z)-Q_{d}(\tau;z)\big{)}\leadsto
N\left(0,\tau(1-\tau)z^{\prime}\left(\mathbb{E}[f_{Y_{d}|X}^{2}(X^{\prime}\theta_{d}(\tau)|X)XX^{\prime}\mathbf{1}\\{D=d\\}]\right)^{-1}z\right).$

What is of interest here is that the asymptotic variance
$\tau(1-\tau)z^{\prime}(\mathbb{E}[f_{Y_{d}|X}^{2}(X^{\prime}\theta_{d}(\tau)|X)XX^{\prime}\mathbf{1}\\{D=d\\}])^{-1}z$
is known to be the semi-parametric efficiency bound for all estimators of the
linear conditional quantile function (Newey and Powell,, 1990). In particular,
the rank-score balanced estimator of the CQF is as efficient as the estimate
of the CQF based on the weighted quantile regression program (Koenker,, 2005;
Koenker and Zhao,, 1994; Zhao,, 2001). This lends further support to the
heuristic arguments made in Section 3.2. Though we note that as the
conditional densities can be hard to estimate, the weighted quantile
regression problems can be less popular in practice.

Since the rank-score debiased estimates of $Q_{1}(\tau;z)$ and $Q_{0}(\tau;z)$
are asympototically independent, Theorem 1 and the Continuous Mapping Theorem
yield the following result for the HQTE process.

######  Theorem 2 (Weak convergence of the rank-score debiased HQTE process).

Let $\mathcal{T}$ be a compact subset of $(0,1)$. Under the conditions of
Theorem 1,

$\displaystyle\sqrt{n}\big{(}\widehat{\alpha}(\cdot;z)-\alpha(\cdot;z)\big{)}\leadsto\mathbb{G}_{1}(\cdot\>;z)+\mathbb{G}_{0}(\cdot\>;z)\quad{}\mathrm{in}\quad{}\ell^{\infty}(\mathcal{T}),$

where $\mathbb{G}_{1}(\cdot\>;z)$, $\mathbb{G}_{0}(\cdot\>;z)$ are
independent, centered Gaussian processes with covariance functions
$(\tau_{1},\tau_{2})\mapsto H_{d}(\tau_{1},\tau_{2};z)$ with $d\in\\{0,1\\}$.

The takeaway from Theorem 2 is that the HQTE process converges weakly to the
sum of two independent centered Gaussian processes. We illustrate Theorem 2
with four examples; for more elaborate applications of process weak
convergence in the context of quantile regression we refer to Belloni et al.,
2019a ; Chao et al., (2017); Angrist et al., (2006); Chernozhukov and
Fernández-Val, (2005).

######  Example 5 (Asymptotic normality of the HQTE estimator).

For fixed quantile $\tau\in\mathcal{T}$, Theorem 2 implies that
$\sqrt{n}\big{(}\widehat{\alpha}(\tau;z)-\alpha(\tau;z)\big{)}$ is
asymptotically normal with mean zero and variance
$\sigma^{2}(\tau;z):=\lim_{n\rightarrow\infty}\sigma_{(n)}^{2}(\tau;z)$, where

$\displaystyle\sigma_{(n)}^{2}(\tau;z)$
$\displaystyle:=\tau(1-\tau)z^{\prime}\left[\Big{(}\pi_{1}\mathbb{E}[f_{Y_{1}|X}^{2}(X^{\prime}\theta_{1}(\tau)|X)XX^{\prime}\mid
D=1]\Big{)}^{-1}\right.$ $\displaystyle\left.\hskip
60.0pt\quad{}+\Big{(}\pi_{0}\mathbb{E}[f_{Y_{0}|X}^{2}(X^{\prime}\theta_{0}(\tau)|X)XX^{\prime}\mid
D=0]\Big{)}^{-1}\right]z,$

where $0<\pi_{1}=1-\pi_{0}=\mathbb{P}\left\\{D=1\right\\}<1$.

######  Example 6 (Joint asymptotic normality of the HQTE estimator at
finitely many quantiles).

Consider a finite collection of quantile levels
$\\{\tau_{1},\ldots,\tau_{K}\\}\subset\mathcal{T}$. Theorem 2 implies that the
collection
$\sqrt{n}\big{(}\widehat{\alpha}(\tau_{j};z)-\alpha(\tau_{j};z)\big{)}$,
$j=1,\ldots,K$, is jointly asymptotically normal with mean zero and covariance
matrix
$\Sigma=(H_{1}(\tau_{j},\tau_{k};z)+H_{0}(\tau_{j},\tau_{k};z))_{j,k=1}^{K}$.

######  Example 7 (Uniform confidence bands for the HQTE curve).

Define
$K(z):=\sup_{\tau\in\mathcal{T}}|\mathbb{G}_{1}(\tau;z)/\sigma(\tau;z)+\mathbb{G}_{0}(\tau;z)/\sigma(\tau;z)|$,
where $\sigma^{2}(\tau;z)$ is the variance from Example 5. Let
$\hat{\kappa}(\alpha;z)$ and $\widehat{\sigma}_{n}^{2}(\tau;z)$ be (uniformly)
consistent estimates of the $\alpha$ quantile of $K(z)$ and
$\sigma^{2}(\tau;z)$, respectively. Then,

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\left\\{\alpha(\tau;z)\in\left[\hat{\alpha}(\tau;z)\pm\hat{\kappa}(\alpha;z)\frac{\widehat{\sigma}_{n}(\tau;z)}{\sqrt{n}}\right],\>\tau\in\mathcal{T}\right\\}=\alpha.$

A consistent estimate $\hat{\kappa}(\alpha;z)$ can be obtained via simulation
based bootstrap, i.e. sampling from
$\widehat{K}(z)=\sup_{\tau\in\mathcal{T}}|\widetilde{\mathbb{G}}_{1}(\tau;z)+\widetilde{\mathbb{G}}_{0}(\tau;z)|$,
where $\widetilde{\mathbb{G}}_{1}(\tau;z)$ and
$\widetilde{\mathbb{G}}_{0}(\tau;z)$ are independent centered Gaussian
processes with covariance functions based on uniformly consistent plug-in
estimates of the operators $(\tau_{1},\tau_{2})\mapsto
H_{d}(\tau_{1},\tau_{2};z)/(\sigma(\tau_{1};z)\sigma(\tau_{2};z))$,
$d\in\\{0,1\\}$.

######  Example 8 (Asymptotic theory for the integrated HQTE curve).

Assessing the HQTE on a specific quantile is often less relevant than
assessing the average HQTE over a continuum of quantile levels $\mathcal{T}$
(e.g., lower, middle, or upper quantiles). In such cases, it is natural to
consider the integrated HQTE. Theorem 2 and the continuous mapping theorem
imply that
$\sqrt{n}\int_{\mathcal{T}}\big{(}\widehat{\alpha}(\tau;z)-\alpha(\tau;z)\big{)}d\tau\leadsto
I(z)$, where
$I(z):=\int_{\mathcal{T}}\mathbb{G}_{1}(\tau,z)d\tau+\int_{\mathcal{T}}\mathbb{G}_{0}(\tau;z)d\tau.$
While the random variable $I(z)$ is not distribution-free, its distribution
can be approximated via re-sampling techniques (Chernozhukov and Fernández-
Val,, 2005).

###  4.4 Duality theory for the rank-score debiasing program

In this section we introduce the dual to the rank-score debiasing program (8)
and explain its pivotal role in the proofs of the weak convergence results in
Sections 4.3. The dual program is also important for constructing uniformly
consistent estimates of the covariance function in Sections 4.5.

Observe that the solution to the rank-score debiasing program (8) can be
written as
$\widehat{w}(\tau;z)=\widehat{w}_{0}(\tau;z)+\widehat{w}_{1}(\tau;z)$, with
the $\widehat{w}_{d}(\tau;z)$’s being the solutions to two independent
optimization problems:

$\displaystyle\widehat{w}_{d}(\tau;z)\in\arg\min_{w\in\mathbb{R}^{n}}\left\\{\sum_{i=1}^{n}w_{i}^{2}\hat{f}_{i}^{-2}(\tau):\>\left\|z-\frac{1}{\sqrt{n}}\sum_{i:D_{i}=d}w_{i}X_{i}\right\|_{\infty}\leq\frac{\gamma_{d}}{n}\right\\},\quad{}\quad{}d\in\\{0,1\\}.$
(12)

These two optimization problems have the following two duals:

$\displaystyle\hat{v}_{d}(\tau;z)\in\arg\min_{v\in\mathbb{R}^{p}}\>\>\>$
$\displaystyle\left\\{\frac{1}{4n}\sum_{i:D_{i}=d}\hat{f}_{i}^{2}(\tau)(X_{i}^{\prime}v)^{2}+z^{\prime}v+\frac{\gamma_{d}}{n}\|v\|_{1}\right\\},\quad{}\quad{}d\in\\{0,1\\}.$
(13)

Provided that strong duality holds, we can estimate the rank-score debiasing
weights $\widehat{w}_{d}(\tau;z)$ by either solving the primal problems (12)
or by solving the dual problems (13) and exploiting the explicit relationship
between primal and dual solutions. To be precise, we have the following
result:

######  Lemma 5 (Dual characterization of the rank-score debiasing program).

  * (i)

Programs (12) and (13) form a primal-dual pair.

  * (ii)

Let $\delta\in(0,1)$ and $d\in\\{0,1\\}$. Suppose that Conditions 3, 5 (i),
and 7 hold. There exists an absolute constant $c_{1}>1$ such that for all
$\gamma_{d}>0$ that satisfy $\gamma_{d}\geq
c_{1}\varphi_{\max}\kappa_{2}^{-1}(\infty)\bar{f}^{2}\|z\|_{2}\sqrt{n\log(p/\delta)}$,
we have with probability at least $1-\delta$, for all $1\leq i\leq n$ and
$\tau\in\mathcal{T}$,

$\displaystyle\widehat{w}_{d,i}(\tau;z)=\begin{cases}-\frac{\hat{f}_{i}^{2}(\tau)}{2\sqrt{n}}X_{i}^{\prime}\hat{v}_{d}(\tau;z),&i\in\\{j:D_{j}=d\\}\\\
0,&i\notin\\{j:D_{j}=d\\},\end{cases}$

where $\widehat{w}_{d}(\tau;z)$ and $\hat{v}_{d}(\tau;z)$ are the solutions to
the programs (12) and (13), respectively.

The important takeaway from Lemma 5 is that, with high probability, for
$\gamma_{d}>0$ large enough, the rank-score balanced estimator (9) has the
following equivalent dual formulation:

$\displaystyle\widehat{Q}_{d}(\tau;z)=z^{\prime}\hat{\theta}_{d}(\tau)-\frac{1}{2n}\sum_{i:D_{i}=d}\hat{f}_{i}(\tau)\big{(}\tau-\mathbf{1}\\{Y_{i}\leq
X_{i}^{\prime}\hat{\theta}_{d}(\tau)\\}\big{)}X_{i}^{\prime}\hat{v}_{d}(\tau;z).$
(14)

Thus, while the original formulation of the rank-score balanced estimator
involves a complicated sum over the rank-score debiasing weights
$\widehat{w}_{1}(\tau;z),\ldots,\widehat{w}_{n}(\tau;z)$, the dual formulation
is a simple linear function of the dual solution
$\hat{v}_{d}(\tau;z)\in\mathbb{R}^{p}$. Therefore, we can expect that (at
least for fixed $\tau$ and $p$) the rank-score debiased estimator can be
approximated by a sum of $n$ independent and identically distributed random
variables. The following non-asymptotic Bahadur-type representation is a
significantly refined version of this statement (holding uniformly in
$\tau\in\mathcal{T}$ and for $p\geq n$). It is key to the weak convergence
results in Section 4.3.

######  Lemma 6 (Bahadur-type representation).

Let $\mathcal{T}$ be a compact subset of $(0,1)$ and $\delta\in(0,1)$. Suppose
that Conditions 1–10 hold with
$\varrho_{n}=\sqrt{(s_{v}+s_{\theta})\log(np/\delta)/n}$ and
$\epsilon_{n}^{2}=O(\sqrt{n}h^{-1}\varrho_{n}^{2}+h^{2})$. In addition,
suppose that $(s_{v}+s_{\theta})^{2}\log^{2}(np/\delta)\log^{2}(n)=o(nh^{2})$,
$h^{2}s_{v}=o(1)$, and $\|z\|_{2}=O(1)$. If
$\lambda_{d}\asymp\sqrt{\varphi_{\max}}\sqrt{n\log(np/\delta)}$ and
$\gamma_{d}\asymp\|z\|_{2}\left(h^{-1}s_{\theta}\log(np/\delta)+\sqrt{n}h^{2}\right)\sqrt{n}$,
then

$\displaystyle\widehat{Q}_{d}(\tau;z)-Q_{d}(\tau;z)$
$\displaystyle\quad{}\quad{}\quad{}=-\frac{1}{2n}\sum_{i:D_{i}=d}f_{Y_{d}|X}(X_{i}^{\prime}\theta_{d}(\tau)|X_{i})\big{(}\tau-\mathbf{1}\\{Y_{i}\leq
X_{i}^{\prime}\theta_{d}(\tau)\\}\big{)}X_{i}^{\prime}v_{d}(\tau;z)+e_{d}(\tau;z),$

where
$v_{d}(\tau;z)=-2\big{(}\mathbb{E}[f_{Y_{d}|X}^{2}\big{(}X^{\prime}\theta_{d}(\tau)|X\big{)}XX^{\prime}\mathbf{1}\\{D=d\\}]\big{)}^{-1}z$,
and, with probability at least $1-\delta$,

$\displaystyle\sup_{\tau\in\mathcal{T}}|e_{n}(\tau;z)|$ $\displaystyle\lesssim
c_{2}\left(\varrho_{n}^{3/2}(\log
n)^{3/4}+h^{2}\varrho_{n}+h^{-1}\varrho_{n}^{2}\right),$

where $c_{2}>0$ depends on
$\bar{f},\underline{f},L_{f},L_{\theta},C_{Q},\kappa_{1}(2),\kappa_{2}(\infty),\varphi_{\max},\|z\|_{2}$.

The upper bound (or: rate) on the remainder term $e_{n}(\tau;z)$ comprises a
parametric and a non-parametric part. The parametric part is
$\varrho_{n}^{3/2}(\log
n)^{3/4}=(s_{v}+s_{\theta})^{3/4}\log^{3/4}(np/\delta)\log^{3/4}(n)/n^{3/4}$.
Up to the $\log$-factors this rate matches the optimal rate of the residuals
of the Bahadur representation for classical estimators of the quantile
function (Bahadur,, 1966; Kiefer,, 1967) as well as quantile regression
estimators in low dimensions (Zhou and Portnoy,, 1996). The non-parametric
part $h^{2}\varrho_{n}+h^{-1}\varrho_{n}^{2}$ depends on the bandwidth $h>0$.
The particular dependence of the bandwidth is the result of the two-fold
dependence of the rank-score debiased estimator on the non-parametric density
estimates: a direct dependence via $\hat{f}_{i}(\tau)$ and an indirect
dependence via $\widehat{v}_{d}(\tau;z)$.

###  4.5 Consistent estimates of the covariance function

The weak convergence results and examples from Section 4.3 are only
practically relevant together with an estimator of the asymptotic covariance
function that is uniformly consistent in $\tau_{1},\tau_{2}\in\mathcal{T}$.
Here, we show how to exploit the duality formalism from Section 4.4 to
construct such estimators.

An estimate for the covariance function $(\tau_{1},\tau_{2})\mapsto
H_{d}(\tau_{1},\tau_{2};z)$ is given by

$\displaystyle\widehat{H}_{d}(\tau_{1},\tau_{2};z):=(\tau_{1}\wedge\tau_{2}-\tau_{1}\tau_{2})\hat{v}_{d}^{\prime}(\tau_{1};z)\left(\frac{1}{4n}\sum_{i:D_{i}=d}\hat{f}_{i}(\tau_{1})\hat{f}_{i}(\tau_{2})X_{i}X_{i}^{\prime}\right)\hat{v}_{d}(\tau_{2};z),$

where $\hat{v}_{d}(\tau;z)$ is the solution to the dual program (13) (see
Section 4.4). By the following lemma this estimate is uniformly consistent in
$\tau_{1},\tau_{2}\in\mathcal{T}$.

######  Lemma 7.

Recall the setup of Theorem 1 and let
$\varrho_{n}=\sqrt{(s_{v}+s_{\theta})\log(np)/n}$. The following holds:

$\displaystyle\sup_{\tau_{1},\tau_{2}\in\mathcal{T}}\left|\widehat{H}_{d}(\tau_{1},\tau_{2};z)-H_{d}(\tau_{1},\tau_{2};z)\right|=O_{p}\left(c_{3}\Big{(}\varrho_{n}+\sqrt{n}\sqrt{s_{v}}h^{-1}\varrho_{n}^{2}+\sqrt{s_{v}}h^{2}\Big{)}\right),$

where $c_{3}>0$ depends on
$\bar{f},\underline{f},L_{f},L_{\theta},C_{Q},\kappa_{1}(2),\kappa_{2}(\infty),\varphi_{\max},\|z\|_{2}$.

As a consequence, a uniformly consistent estimate of the asymptotic variance
$\sigma^{2}(\tau;z)$ of the HQTE process at a single quantile
$\tau\in\mathcal{T}$ (see Example 5) is given by

$\displaystyle\widehat{\sigma}_{1}^{2}(\tau;z):=\tau(1-\tau)\left(\frac{1}{4n}\sum_{i:D_{i}=1}\hat{f}_{i}^{2}(\tau)\big{(}X_{i}^{\prime}\hat{v}_{1}(\tau;z)\big{)}^{2}+\frac{1}{4n}\sum_{i:D_{i}=0}\hat{f}_{i}^{2}(\tau)\big{(}X_{i}^{\prime}\hat{v}_{0}(\tau;z)\big{)}^{2}\right),$
(15)

where $\hat{v}_{1}(\tau;z)$ and $\hat{v}_{0}(\tau;z)$ are the solutions to the
dual problems (13). The duality formalism from Section 4.4 implies that
another uniformly consistent estimate for $\sigma^{2}(\tau;z)$ is given by

$\displaystyle\widehat{\sigma}_{2}^{2}(\tau;z):=\tau(1-\tau)\sum_{i=1}^{n}\widehat{w}_{i}^{2}(\tau;z)\hat{f}_{i}^{-2}(\tau),$
(16)

where the $\widehat{w}_{i}(\tau;z)$’s are the rank-score debiasing weights.
Neither of the two estimates requires inverting a (high-dimensional) matrix,
which may be surprising given the form of the target $\sigma^{2}(\tau;z)$.

##  5 A practical guide to the rank-score debiasing procedure

In the following, we explain how we implement the rank-score debiasing
procedure with the help of the dual problem. As the rank-score debiasing
estimator of the HQTE depends on the four regularization parameters
$\lambda_{0}$, $\lambda_{1}$, $\gamma_{0}$, and $\gamma_{1}>0$ and the
bandwidth $h>0$ of the non-parametric density estimator, we also explain how
to choose these parameters in robust and data-dependent ways.

###  5.1 Implementing the $\ell_{1}$-penalized quantile regression program

To select $\lambda_{d}>0$ in a data dependent way, we substantially deviate
from the vanilla quantile regression program (6) and instead implement the
weighted $\ell_{1}$-penalized quantile regression problem by Belloni and
Chernozhukov, (2011). That is, we compute the pilot estimate of
$\theta_{d}(\tau)$ as

$\displaystyle\hat{\theta}_{d}(\tau)\in\arg\min_{\theta\in\mathbb{R}^{p}}\left\\{\sum_{i:D_{i}=d}\rho_{\tau}(Y_{i}-X_{i}^{\prime}\theta)+\lambda_{d}\sqrt{\tau(1-\tau)}\sum_{k=1}^{p}\widehat{\sigma}_{d,k}|\theta_{k}|\right\\},$
(17)

with $\widehat{\sigma}_{d,k}^{2}=n^{-1}\sum_{i;D_{i}=d}X_{ik}^{2}$ and
$\lambda_{d}=1.5\cdot\Lambda_{d}(0.9|X_{1},\ldots,X_{n})$, where
$\Lambda_{d}(0.9|X_{1},\ldots,X_{n})$ is the 90%-quantile of
$\Lambda_{d}|X_{1},\ldots,X_{n}$ and

$\displaystyle\Lambda_{d}:=\sup_{\tau\in\mathcal{T}}\max_{1\leq k\leq
p}\left|\sum_{i:D_{i}=d}\frac{(\tau-1\\{U_{i}\leq\tau\\})X_{ik}}{\widehat{\sigma}_{d,k}\sqrt{\tau(1-\tau)}}\right|,$

with $U_{1},\ldots,U_{n}$ be i.i.d. Uniform(0,1) random variables, independent
of $X_{1},\ldots,X_{n}$

###  5.2 Implementing the rank-score debiasing program

Recall the primal and dual programs (12) and (13), respectively. Provided that
strong duality holds, we can estimate the rank-score balancing weights
$\widehat{w}(\tau;z)$ by solving either of the two problems. However, from a
statistical and computational point of view, it is preferable to solve the
dual problems.

First, since the dual programs (13) are unconstrained optimization problems
they allow us to choose the tuning parameter $\gamma_{d}>0$ systematically via
cross-validation. In contrast, the primal problems are constrained
optimization problems which do not naturally lend themselves to cross-
validation procedures. In the simulation study we therefore implement a
10-fold cross validation procedure on the dual problems and choose
$\gamma_{d}>0$ as the smallest tuning parameter which yields a risk that is at
most one standard deviation away from the smallest cross-validated risk. The
main point of this one-standard-deviation (1SE) rule is to estimate debiasing
weights with small bias
$\|z-\frac{1}{\sqrt{n}}\sum_{i:D_{i}=d}w_{i}X_{i}\|_{\infty}$ whose risk is
comparable to the one of the optimal weights. A smaller $\gamma_{d}>0$
produces a less biased estimate, which leads to a better coverage probability
of the confidence interval. It is instructive to compare our 1SE rule with the
1SE rule popularized by Breiman et al., (1984). Breiman et al., (1984) aim to
improve the out-of-sample (classification) accuracy of their estimator and
hence advocate choosing the least variable model whose risk is comparable to
the model with the smallest cross-validated risk. In contrast, we aim to
improve statistical inferential validity and hence are less concerned about
the variability of our estimate than its bias.

Second, since the primal problems (12) are constrained optimization programs,
finding feasible points can be difficult. In contrast, the dual programs are
unconstrained convex optimization problems and therefore can be easily solved
by off-shelf optimization packages. In our simulation studies, we solve the
primal problem using R package CVXR (Fu et al.,, 2017), and the dual problem
using Alternating Direction Method of Multipliers (ADMM) to the
$l_{1}$-regularized quadratic program (Wahlberg et al.,, 2012) via R package
accSDA (Atkins et al.,, 2017).

Third, since the dual programs do not involve the inverses of the estimated
densities $\hat{f}_{i}(\tau)$ they are numerically more stable than the primal
problems. Therefore, in the simulation studies and the real data analysis we
only report results obtained via the dual problem.

###  5.3 Selecting bandwidth $h$ for the non-parametric density estimator

To stabilize the density estimator (7), we replace the $\ell_{1}$-penalized
quantile regression estimates with refitted quantile regression estimates. The
refitted estimates are obtained by fitting a quantile regression model to the
data using only the covariates in the support set of $\hat{\theta}_{d}(\tau)$.
As this density estimator takes a similar form as the one in Belloni et al.,
2019b , we follow their advice and set bandwidth
$h=\min\\{n^{-1/6},\tau(1-\tau)/2\\}$.

##  6 Simulation study

We carry out simulation studies to investigate the performance of the rank-
score debiased estimator. The goal of the simulation studies is to: (1)
illustrate our rank-score debiased estimator provides consistent estimate of
HQTE with nominal-level coverage probabilities, (2) showcase the rank-score
debiased estimator is more efficient than the unweighted quantile regression
estimator, and (3) provide numerical evidence supporting the theoretical
results from Section 4.3.

###  6.1 Simulation design

Our simulation design mimics high-dimensional observation studies where
treatments are assigned based on covariates. We consider the following
generative model:

$\displaystyle Y_{1}=X^{\prime}\theta_{1}+\varepsilon\sigma_{1}(X),\qquad
Y_{0}=X^{\prime}\theta_{0}+\varepsilon\sigma_{0}(X),\qquad
X\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\varepsilon,\qquad\varepsilon\sim N(0,1),$
$\displaystyle D\mid
X\sim\text{Bernoulli}\left(\frac{e^{1-X_{7}+X_{8}}}{1+e^{1-X_{7}+X_{8}}}\right),\qquad
Y=DY_{1}+(1-D)Y_{0}.$

For the noise level $\sigma_{d}(X)$ and the covariates $X$, we consider two
sets of covariate designs for the homoscedastic case and the heteroscedastic
case. We first generate $W\sim N(0,\Sigma)$ where
$\Sigma=(\Sigma_{jk})_{j,k=1}^{p-1}$ and $\Sigma_{jk}=0.5^{|j-k|}$. Then, in
the homoscedastic case, we set $\sigma_{1}(X)=\sigma_{0}(X)=1$ and generate
the covariates with $X_{1}=1$ and $X_{j}=W_{j}$, for $2\leq j\leq p$. In the
heteroscedastic case, we set $X_{1}=1$, $X_{2}=|W_{2}|+0.1$,
$X_{3}=W_{3}^{2}+0.5$, $X_{j}=W_{j}$ for $4\leq j\leq p$, and
$\sigma_{d}(X)=(1-d)X_{2}+dX_{3}$ for $d\in\\{0,1\\}$. In both cases, we set
$\theta_{0}=(0.5,0,1,-1,0,\ldots,0)^{\prime}\in\mathbb{R}^{p}$ and consider
the following three scenarios for $\theta_{1}$: sparse
($\theta_{1}\propto(1,1,1,1,1,1,0,\ldots,0)^{\prime}$), dense
($\theta_{1}\propto(1,1/\sqrt{2},\ldots,1/\sqrt{p})^{\prime}$) and pseudo
dense ($\theta_{1}\propto(1,1/{2},\ldots,1/{p})^{\prime}$). We consider three
different signal strengths $||\theta_{1}||_{2}\in\\{1,2,4\\}$. We choose the
sample size $n$ and the dimension of the covariates $p$ from
$(n,p)\in\\{(600,400),(1000,600)\\}$. As we estimate the CQF separately by
using the observed data in the treated and control groups, the effective
sample size for our rank-score debiasing program $n_{d}$ is approximately half
of the sample size. Thus, the effective sample size is always less than $p$.
Lastly, we set $z=(0,1/\sqrt{2},1/\sqrt{2},0,\ldots,0)^{\prime}$ or
$z=(1,1,1/\sqrt{2},\ldots,1/\sqrt{p})^{\prime}$. Under this data generating
process, the HQTE at $z$ is the linear function
$\alpha(\tau;z)=z^{\prime}(\theta_{1}(\tau)-\theta_{0}(\tau))$.

We implement the rank-score debiased estimator as discussed in Section 5. In
particular, this means that even in the case of homoscedastic noise we do not
use a specialized density estimator that could exploit this extra information.
Since in practice homoscedasticity may be difficult to detect, we do not to
want rely on the validity of the homoscedasticity assumption. To illustrate
the bias-variance trade-off that underlies the tuning parameter $\gamma_{d}$
we report results not just for the “1SE” rule (“Rank-1SE”) but also for a
“2SE” rule (“Rank-2SE”). The “2SE” rule chooses the smallest $\gamma_{d}>0$
that is less than two standard errors away from the tuning parameter with the
lowest dual loss function.

To showcase the merit of the rank-score debiased estimator, we compare it with
the following four methods: “Unweighted Oracle”, “Refit”, “Lasso” and
“Debiased. The “Unweighted Oracle” method fits a quantile regression model
based on the true model, i.e. based on the covariates in the support set of
$\theta_{d}(\tau)$ only. The (unweighted) oracle estimate of the HQTE is
$\widehat{\alpha}^{\mathrm{oracle}}(\tau;z)=z^{\prime}(\hat{\theta}_{1}^{\mathrm{oracle}}(\tau)-\hat{\theta}_{0}^{\mathrm{oracle}}(\tau))$.
We compute this (unweighted) oracle estimate only in the scenario with sparse
$\theta_{d}(\tau)$. The “Refit” method is the following two-step procedure: We
first obtain estimates $\hat{\theta}_{d}(\tau)$ by solving the weighted
$\ell_{1}$-penalized quantile regression program (17). Then, we compute a
refitted estimate $\hat{\theta}_{d}^{\mathrm{refit}}(\tau)$ by fitting a
quantile regression model based only on the covariates in the support set of
$\hat{\theta}_{d}(\tau)$. The refitted estimate of the HQTE is
$\widehat{\alpha}^{\mathrm{refit}}(\tau;z)=z^{\prime}(\hat{\theta}_{1}^{\mathrm{refit}}(\tau)-\hat{\theta}_{0}^{\mathrm{refit}}(\tau))$.
The “Lasso” method refers to simply using the estimates
$\hat{\theta}_{d}(\tau)$ from the $\ell_{1}$-penalized quantile regression
program (17) without any adjustments. The Lasso estimate of the HQTE is thus
$\widehat{\alpha}^{\mathrm{lasso}}(\tau;z)=z^{\prime}(\hat{\theta}_{1}^{\mathrm{lasso}}(\tau)-\hat{\theta}_{0}^{\mathrm{lasso}}(\tau))$.
The “Debiased” method refers to the debiased $\ell_{1}$-penalized quantile
regression coefficient estimate proposed by Zhao et al., (2019). We denote the
debiased estimate of the HQTE as
$\widehat{\alpha}^{\mathrm{debias}}(\tau;z)=z^{\prime}(\hat{\theta}_{1}^{\mathrm{debias}}(\tau)-\hat{\theta}_{0}^{\mathrm{debias}}(\tau))$
with

$\displaystyle\hat{\theta}_{d}^{\mathrm{debias}}(\tau)=\hat{\theta}_{d}^{\mathrm{Lasso}}(\tau)+\hat{\Theta}_{d}(\tau)\cdot\frac{1}{n_{d}}\sum_{i:D_{i}=d}X_{i}(\tau-\mathbf{1}\\{Y_{i}\leq
X_{i}^{\prime}\hat{\theta}_{d}(\tau)\\}),$

where $\hat{\Theta}_{d}(\tau)$ is an estimate of the inverse covariance matrix
$\big{[}\mathbb{E}[f_{Y_{d}|X}(X^{\prime}\theta_{d}(\tau)|X)XX^{\prime}\mathbf{1}\\{D=d\\}]\big{]}^{-1}$.
Following the recommendation by Zhao et al., (2019), we use the R package
clime to obtain $\hat{\Theta}_{d}(\tau)$.

Confidence intervals for the rank-score debiased estimator are based on the
asymptotic normality results in Section 4.3 and hold under the mild regularity
conditions stated in Section 4.1. Confidence intervals for the (unweighted)
oracle method are constructed using standard large sample theory (Angrist et
al.,, 2006). Confidence intervals for the Lasso and the Refit method are
constructed assuming that the selected models equal the true model, i.e. the
support sets of $\hat{\theta}_{d}^{\mathrm{lasso}}(\tau)$,
$\hat{\theta}_{d}^{\mathrm{refit}}(\tau)$ equal the support set of
$\theta_{d}(\tau)$. This assumption is satisfied under strong oracle
conditions (Fan and Li,, 2001). As Zhao et al., (2019) focus on providing
accurate point estimate of the quantile regression coefficients, they do not
construct confidence interval for $\theta_{d}(\tau)$. Based on our conjecture
provided in the Supplementary Materials, we construct confidence intervals
based on normal approximation with an estimated asymptotic variance equal to
$\pi_{1}z^{\prime}\hat{\Theta}_{1}(\tau)\frac{1}{n_{1}}\sum_{i:D_{i}=1}X_{i}X_{i}^{\prime}\hat{\Theta}_{1}(\tau)z+\pi_{0}z^{\prime}\hat{\Theta}_{0}(\tau)\frac{1}{n_{0}}\sum_{i:D_{i}=0}X_{i}X_{i}^{\prime}\hat{\Theta}_{0}(\tau)z$.

Figure 1: Simulation results for homoscedastic data with $(n,p)=(1000,600)$,
sparse $\theta_{1}$ with $||\theta_{1}||_{2}=4$ and sparse $z$. Panel (A) Bias
comparison. Panel (B) Variance comparison. Panel (C) Coverage probability of
the confidence interval while the nominal coverage probability is 95%. Panel
(D) Histograms of the standardized estimates of the rank-score debiased
estimator displayed in (18), density of $N(0,1)$ in red. Figure 2: Simulation
results for heteroscedastic data with $(n,p)=(1000,600)$, sparse $\theta_{1}$
with $||\theta_{1}||_{2}=4$, and sparse $z$. Panel (A) Bias comparison. Panel
(B) Variance comparison. Panel (C) Coverage probability of the confidence
interval while the nominal coverage probability is 95%. Panel (D) Histograms
of the standardized estimates of the rank-score debiased estimator displayed
in (18), density of $N(0,1)$ in red.

###  6.2 Simulation results

We measure the performance of the estimators in terms of their biases
(computed as the differences between the mean of the Monte Carlo estimates of
$\alpha(\tau;z)$ and the true HQTE), variances (computed as the variances of
the Monte Carlo estimates of $\alpha(\tau;z)$) and coverage probabilities of
the confidence intervals with the nominal coverage probability of 95%. We
provide finite sample comparisons through Table 1 and Figure 1 for
homoscedastic data, and Table 2 and Figure 2 for heteroscedastic data. Details
about the model parameters are given in the captions of these tables and
figures. Our simulation results are evaluated through 2,000 Monte Carlo
samples.

The main takeaway from the simulation study is that the rank-score debiased
estimator with $\gamma_{d}$ selected by the 1SE rule outperforms the Refitted
and Lasso estimators in terms of bias, variance and the validity of inference
in most scenarios. In the following we highlight three conclusions. First, the
rank-score debiased estimator performs better in sparse than in dense models.
Second, the rank-score debiased estimator can have a smaller variance than the
unweighted Oracle estimator in the heteroscedastic cases when $\theta_{1}$ is
sparse. Third, the asymptotic normality results from Section 4.3 continue to
hold reasonably well in finite samples. This can be deduced from Figures 1(D)
and 2(D), in which we provide histograms of the standardized estimates of the
rank-score debiased estimator,

$\displaystyle\widehat{\sigma}_{2}^{-1}(\tau;z)\cdot\sqrt{n}\big{(}\widehat{\alpha}(\tau;z)-\alpha(\tau;z)\big{)},$
(18)

where $\widehat{\sigma}_{2}(\tau;z)$ is the estimate defined in eq. (15).
These histograms fit the overlaid N(0,1) densities. In contrast, the Lasso
estimator is clearly biased (Figures 1(A) and 2(B)) and so is the Refitting
estimator in scenarios with smaller signal to noise ratio, i.e. scenarios with
small $\|\theta_{1}(\tau)\|_{2}$ (Tables 1 and 2). These biases suggest that
the oracle condition is violated and hence the finite sample distributions of
these estimators may not be approximated by a standard normal distribution.
The debiased quantile Lasso estimator has small biases, but it often has
larger variance compared to the rank-score debiased estimator. This
observation is in-line with our conjecture based on the derivation provided in
the Supplementary Materials.

| $\tau$ | Unweighted Oracle | Refit | Lasso | Debias | Rank-1SE | Rank-2SE
---|---|---|---|---|---|---|---
|  |  $n=600,p=400$, sparse $\theta_{1}$ with $||\theta_{1}||_{2}=1$, sparse $z$
$\sqrt{n}$Bias | $0.2$ | $0.12(0.08)$ | $-0.03(0.08)$ | $-1.25(0.08)$ | $-0.37(0.10)$ | $-0.50(0.09)$ | $-0.47(0.09)$
| $0.5$ | $0.14(0.07)$ | $-0.1(0.07)$ | $-1.15(0.08)$ | $-0.41(0.09)$ | $-0.27(0.08)$ | $-0.23(0.09)$
| $0.7$ | $0.10(0.07)$ | $-0.11(0.07)$ | $-1.16(0.07)$ | $-0.54(0.08)$ | $-0.26(0.08)$ | $-0.25(0.08)$
$n$Var | $0.2$ | $6.58(0.31)$ | $7.10(0.33)$ | $6.48(0.35)$ | $9.58(0.47)$ | $7.97(0.37)$ | $8.06(0.38)$
| $0.5$ | $5.13(0.21)$ | $5.60(0.24)$ | $6.13(0.34)$ | $7.37(0.35)$ | $6.53(0.33)$ | $6.77(0.33)$
| $0.7$ | $5.14(0.21)$ | $5.56(0.25)$ | $5.13(0.28)$ | $7.14(0.29)$ | $5.9(0.25)$ | $6.01(0.26)$
Coverage | $0.2$ | $0.94$ | $0.85$ | $0.84$ | $0.97$ | $0.91$ | $0.93$
| $0.5$ | $0.95$ | $0.91$ | $0.85$ | $0.98$ | $0.95$ | $0.96$
| $0.7$ | $0.95$ | $0.89$ | $0.88$ | $0.96$ | $0.95$ | $0.95$
|  |  $n=600,p=400$, pseudo sparse $\theta_{1}$ with $||\theta_{1}||_{2}=1$, sparse $z$
$\sqrt{n}$Bias | $0.2$ | $-$ | $-0.34(0.09)$ | $-2.15(0.08)$ | $-1.20(0.10)$ | $-0.90(0.09)$ | $-0.86(0.09)$
| $0.5$ | $-$ | $-0.33(0.08)$ | $-1.60(0.08)$ | $-0.93(0.09)$ | $-0.64(0.09)$ | $-0.58(0.09)$
| $0.7$ | $-$ | $-0.27(0.08)$ | $-1.54(0.08)$ | $-1.09(0.09)$ | $-0.75(0.08)$ | $-0.72(0.09)$
$n$Var | $0.2$ | $-$ | $7.55(0.32)$ | $6.75(0.51)$ | $9.15(0.52)$ | $8.03(0.41)$ | $8.05(0.4)$
| $0.5$ | $-$ | $6.39(0.27)$ | $6.46(0.41)$ | $8.61(0.49)$ | $8.36(0.43)$ | $8.96(0.46)$
| $0.7$ | $-$ | $6.15(0.31)$ | $5.85(0.37)$ | $8.21(0.42)$ | $7.01(0.32)$ | $7.46(0.34)$
Coverage | $0.2$ | $-$ | $0.82$ | $0.73$ | $0.97$ | $0.97$ | $0.98$
| $0.5$ | $-$ | $0.88$ | $0.82$ | $0.96$ | $0.94$ | $0.96$
| $0.7$ | $-$ | $0.86$ | $0.77$ | $0.97$ | $0.96$ | $0.98$
|  |  $n=600,p=400$, dense $\theta_{1}$ with $||\theta_{1}||_{2}=1$, sparse $z$
$\sqrt{n}$Bias | $0.2$ | $-$ | $-2.02(0.12)$ | $-3.06(0.09)$ | $-2.15(0.12)$ | $-1.91(0.10)$ | $-1.79(0.1)$
| $0.5$ | $-$ | $-1.76(0.11)$ | $-3.04(0.09)$ | $-1.73(0.12)$ | $-1.55(0.10)$ | $-1.46(0.10)$
| $0.7$ | $-$ | $-1.64(0.11)$ | $-2.94(0.09)$ | $-2.14(0.12)$ | $-1.83(0.09)$ | $-1.73(0.10)$
$n$Var | $0.2$ | $-$ | $15.30(0.91)$ | $7.39(0.74)$ | $13.80(0.83)$ | $9.72(0.60)$ | $9.96(0.59)$
| $0.5$ | $-$ | $12.55(0.75)$ | $8.45(0.75)$ | $14.6(0.79)$ | $9.99(0.54)$ | $10.38(0.54)$
| $0.7$ | $-$ | $12.45(0.77)$ | $7.59(0.70)$ | $14.09(0.83)$ | $9.06(0.55)$ | $9.14(0.53)$
Coverage | $0.2$ | $-$ | $0.65$ | $0.63$ | $0.96$ | $0.93$ | $0.94$
| $0.5$ | $-$ | $0.70$ | $0.60$ | $0.98$ | $0.96$ | $0.98$
| $0.7$ | $-$ | $0.72$ | $0.60$ | $0.96$ | $0.97$ | $0.98$
|  |  $n=600,p=400$, sparse $\theta_{1}$ with $||\theta_{1}||_{2}=2$, dense $z$
$\sqrt{n}$Bias | $0.2$ | $0.79(0.12)$ | $0.77(0.13)$ | $4.79(0.11)$ | $2.00(0.15)$ | $2.71(0.12)$ | $2.64(0.12)$
| $0.5$ | $0.44(0.11)$ | $0.35(0.11)$ | $-1.32(0.12)$ | $-0.81(0.15)$ | $-0.26(0.13)$ | $-0.24(0.13)$
| $0.7$ | $0.57(0.11)$ | $0.56(0.13)$ | $-1.16(0.11)$ | $-1.03(0.16)$ | $-0.49(0.13)$ | $-0.49(0.13)$
$n$Var | $0.2$ | $14.59(0.69)$ | $16.52(0.73)$ | $11.33(1.49)$ | $23.32(1.27)$ | $14.62(0.99)$ | $14.89(0.98)$
| $0.5$ | $11.22(0.52)$ | $12.55(0.58)$ | $14.97(0.73)$ | $22.43(1.09)$ | $16.69(0.82)$ | $16.97(0.84)$
| $0.7$ | $11.66(0.55)$ | $16.95(0.83)$ | $12.61(0.62)$ | $27.38(1.27)$ | $17.63(0.87)$ | $17.97(0.89)$
Coverage | $0.2$ | $0.95$ | $0.95$ | $0.78$ | $0.99$ | $0.93$ | $0.94$
| $0.5$ | $0.95$ | $0.93$ | $0.88$ | $0.99$ | $0.95$ | $0.97$
| $0.7$ | $0.97$ | $0.88$ | $0.94$ | $0.98$ | $0.96$ | $0.96$

Table 1: Homoscedastic data. Standard errors of estimates based on 1,000 Monte
Carlo samples are given in parenthesis. All standard errors of the coverage
probability are smaller than 0.01 and thus are omitted.

| $\tau$ | Unweighted Oracle | Refit | Lasso | Debias | Rank-1SE | Rank-2SE
---|---|---|---|---|---|---|---
|  |  $n=600,p=400$, sparse $\theta_{1}$ with $||\theta_{1}||_{2}=1$, sparse $z$
$\sqrt{n}$Bias | $0.2$ | $0.21(0.18)$ | $-10.65(0.13)$ | $-7.77(0.07)$ | $-1.14(0.15)$ | $-3.71(0.1)$ | $-3.51(0.1)$
| $0.5$ | $-0.29(0.12)$ | $-2.1(0.15)$ | $-5.76(0.11)$ | $-0.97(0.16)$ | $-1.71(0.14)$ | $-1.22(0.15)$
| $0.7$ | $0.2(0.15)$ | $0.24(0.16)$ | $-2.59(0.14)$ | $-1.22(0.15)$ | $-0.83(0.13)$ | $-0.73(0.13)$
$n$Var | $0.2$ | $31.10(1.33)$ | $16.39(4.63)$ | $4.42(2.26)$ | $24.18(1.16)$ | $10.26(0.98)$ | $10.81(0.95)$
| $0.5$ | $15.54(0.61)$ | $23.07(1.02)$ | $11.26(1.82)$ | $24.29(1.10)$ | $17.34(0.91)$ | $17.80(0.92)$
| $0.7$ | $22.25(0.85)$ | $26.19(1.09)$ | $21.13(1.32)$ | $23.79(1.09)$ | $17.66(0.81)$ | $18.01(0.82)$
Coverage | $0.2$ | $0.92$ | $0.17$ | $0.12$ | $0.97$ | $0.88$ | $0.90$
| $0.5$ | $0.94$ | $0.64$ | $0.50$ | $0.96$ | $0.90$ | $0.92$
| $0.7$ | $0.94$ | $0.91$ | $0.90$ | $0.99$ | $0.94$ | $0.98$
|  |  $n=600,p=400$, pseudo sparse $\theta_{1}$ with $||\theta_{1}||_{2}=1$, sparse $z$
$\sqrt{n}$Bias | $0.2$ | $-$ | $-7.75(0.13)$ | $-4.55(0.08)$ | $1.00(0.17)$ | $-0.90(0.11)$ | $-0.46(0.12)$
| $0.5$ | $-$ | $-2.79(0.14)$ | $-4.69(0.09)$ | $-1.44(0.14)$ | $-1.81(0.12)$ | $-1.38(0.13)$
| $0.7$ | $-$ | $-0.63(0.17)$ | $-3.00(0.16)$ | $-2.92(0.17)$ | $-1.76(0.15)$ | $-1.50(0.16)$
$n$Var | $0.2$ | $-$ | $16.79(2.82)$ | $5.84(1.02)$ | $29.20(1.49)$ | $14.53(0.49)$ | $17.20(0.59)$
| $0.5$ | $-$ | $18.50(1.16)$ | $8.92(1.31)$ | $19.34(1.10)$ | $13.91(0.88)$ | $15.44(0.95)$
| $0.7$ | $-$ | $30.67(1.43)$ | $25.12(1.55)$ | $29.61(1.78)$ | $23.08(1.28)$ | $24.71(1.25)$
Coverage | $0.2$ | $-$ | $0.26$ | $0.39$ | $0.96$ | $0.96$ | $0.97$
| $0.5$ | $-$ | $0.63$ | $0.54$ | $0.96$ | $0.93$ | $0.94$
| $0.7$ | $-$ | $0.84$ | $0.85$ | $0.98$ | $0.95$ | $0.96$
|  |  $n=600,p=400$, dense $\theta_{1}$ with $||\theta_{1}||_{2}=1$, sparse $z$
$\sqrt{n}$Bias | $0.2$ | $-$ | $-10.94(0.16)$ | $-3.59(0.07)$ | $-0.28(0.19)$ | $-1.38(0.13)$ | $-1.31(0.14)$
| $0.5$ | $-$ | $-3.53(0.15)$ | $-4.60(0.10)$ | $-1.10(0.17)$ | $-1.49(0.14)$ | $-1.43(0.14)$
| $0.7$ | $-$ | $-3.61(0.20)$ | $-4.39(0.16)$ | $-4.41(0.21)$ | $-4.02(0.18)$ | $-3.93(0.18)$
$n$Var | $0.2$ | $-$ | $25.23(6.01)$ | $4.87(0.77)$ | $35.73(1.44)$ | $18.35(0.85)$ | $19.27(0.86)$
| $0.5$ | $-$ | $24.13(1.45)$ | $10.84(1.24)$ | $29.81(1.60)$ | $19.78(0.95)$ | $20.12(0.95)$
| $0.7$ | $-$ | $41.7(2.54)$ | $27.39(2.00)$ | $43.87(2.91)$ | $31.38(2.24)$ | $31.77(2.23)$
Coverage | $0.2$ | $-$ | $0.03$ | $0.61$ | $0.96$ | $0.93$ | $0.94$
| $0.5$ | $-$ | $0.60$ | $0.57$ | $0.99$ | $0.95$ | $0.97$
| $0.7$ | $-$ | $0.72$ | $0.72$ | $0.98$ | $0.93$ | $0.94$
|  |  $n=600,p=400$, sparse $\theta_{1}$ with $||\theta_{1}||_{2}=2$, dense $z$
$\sqrt{n}$Bias | $0.2$ | $0.32(0.10)$ | $-6.47(0.12)$ | $-8.68(0.20)$ | $-7.63(0.19)$ | $-6.37(0.19)$ | $-6.23(0.19)$
| $0.5$ | $0.36(0.07)$ | $0.08(0.08)$ | $-1.71(0.08)$ | $2.08(0.10)$ | $-0.55(0.08)$ | $-0.46(0.08)$
| $0.7$ | $0.64(0.1)$ | $0.96(0.11)$ | $1.78(0.10)$ | $2.82(0.12)$ | $0.43(0.10)$ | $0.47(0.10)$
$n$Var | $0.2$ | $13.06(0.58)$ | $15.93(1.03)$ | $38.42(4.95)$ | $38.25(4.39)$ | $34.45(4.04)$ | $34.39(3.97)$
| $0.5$ | $5.85(0.22)$ | $5.74(0.28)$ | $6.69(0.46)$ | $9.11(0.64)$ | $6.27(0.33)$ | $6.43(0.33)$
| $0.7$ | $11.23(0.48)$ | $13.24(0.64)$ | $10.16(0.62)$ | $13.7(1.00)$ | $10.71(0.51)$ | $10.78(0.52)$
Coverage | $0.2$ | $0.96$ | $0.60$ | $0.39$ | $0.66$ | $0.82$ | $0.85$
| $0.5$ | $0.96$ | $0.95$ | $0.86$ | $0.99$ | $0.94$ | $0.94$
| $0.7$ | $0.95$ | $0.90$ | $0.92$ | $0.94$ | $0.95$ | $0.95$

Table 2: Heteroscedastic data. Standard errors of estimates based on 1,000
Monte Carlo samples are given in parenthesis. All standard errors of the
coverage probability are smaller than 0.01 and thus are omitted.

##  7 A case study

###  7.1 Study design

To illustrate the advantages of considering HQTE, we apply the proposed method
to study the heterogeneous effect of statin usage, especially when combined
with a healthy lifestyle, in lowering the low–density lipoprotein (LDL)
chlestrol concentration levels for older Alzheimer’s disease (AD) patients
enrolled in the UK Biobank study.

Alzheimer’s disease (AD) is the sixth leading cause of death in the United
States, directly affecting an estimated 5.8 million Americans and incurring
nearly $236 billion of total healthcare costs (Alzheimer’s Association,,
2019). While there is no disease-modifying treatment available for AD, several
studies have reported a reduced risk for progression of AD in statin-treated
populations (Jick et al.,, 2000; Rockwood et al.,, 2002; Geifman et al.,,
2017)111Statin is a commonly prescribed drug due to its clear benefits in
reducing the level of LDL cholesterol through 3-hydroxy-3-methylglutaryl-
coenzyme A reductase (HMGCR) inhibition (Nissen et al.,, 2005).. This slowed
progression of AD might be linked to the reduced cholesterol generation after
statin usage suggested by a substantial body of cellular and molecular
mechanistic evidence (Di Paolo and Kim,, 2011; McGuinness et al.,, 2010).
Thus, our study may provide some additional evidence for the conjecture that
statins are helpful for AD patients as they lower the LDL cholesterol levels.

On the top of management of diseases related to high LDL cholesterol
concentrations, there has been increased global attention on prevention and
risk reduction of AD by maintenance of healthy lifestyle patterns (Barthold et
al.,, 2020; WHO,, 2019; Lourida et al.,, 2019). In this case study, we first
adopt our proposed method to examine if the combined effect of healthy dietary
patterns, increased physical activities, reduced alcohol intake, and reduced
smoking on lowering LDL cholesterol concentration in the statin-treated group
is different from the statin-controlled group in AD patients. Since that the
Mediterranean diet is one of the dietary patterns most commonly investigated
(Kivipelto et al.,, 2018), we define the healthy dietary pattern on the basis
of adherence to the following characteristics: consumption of an increased
amount of fruits, vegetables, fish, and a reduced amount of processed meats
and unprocessed red meats. We then move on to examine if the statin usage
takes heterogeneous effects across different individuals with different
lifestyles in the study cohorts. Detailed descriptions of our data structure
and scientific questions are provided in the next section.

Since existing clinical studies suggest that investigating the benefit of
statin usage can be susceptible to unmeasured confounding factors which induce
potential selection bias, we adopt a genetic variant rs12916-T as a surrogate
treatment variable. This means that if the subject carries the variant
rs12916-T, the treatment indicator variable is set to be one $D=1$, otherwise
is set to be zero. We adopt this genetic surrogate biomarker as the treatment
because the rs12916-T allele only affects the LDL cholesterol concentration
through HMGCR inhibition and is functionally equivalent to statin usage
(Swerdlow et al.,, 2015). Moreover, given that genetic variants are randomly
inherited from parents and are fixed at conception, our treatment variable
(whether the individual carrying rs12916-T) is thus independent of unmeasured
factors such as lifestyle modification after statin usage (Swerdlow et al.,,
2015; Würtz et al.,, 2016). To further avoid potential confounding issues
introduced by genetic pleiotropy and linkage disequilibrium, our approach
includes $p=637$ additional SNPs and lifestyle factors associated with LDL
cholesterol concentration as covariates. We hope that such a study design
makes Condition 1 more plausible to believe in this case study. See
Supplementary Materials for our data pre-processing steps.

Lastly, we recognize that our study design has some potential limitations.
Since the treatment variable is defined as carrying rs12916-T allele or not,
it is a surrogate measurement of statin usage. This suggests that generalizing
the current study findings still warrants further confirmation from clinical
trials.

Figure 3: Panel (A) LDL plasma concentration for the treated and control
groups. Panel (B) HQTE of statin usage by healthy lifestyles. Uniform 95%
confidence bands discussed in Example 7 are given by the shaded regions.

###  7.2 Data structure and analysis

Our data source is the UK Biobank study. The UK Biobank study cohort is a
prospective cohort that enrolled about 500,000 individuals aged from 40 to 69
in the United Kingdom, started in 2006. We focus on AD (and AD proxy) patients
older than 65 years since the majority of AD patients experience their first
symptoms in their mid-60 (Jack Jr et al.,, 2010). To avoid complications due
to missing data, we only include patients with complete covariates. This
results in a cohort of $n_{0}=\sum_{i=1}^{n}(1-D_{i})=563$ subjects who do not
carry the variant rs12916-T, and $n_{1}=\sum_{i=1}^{n}D_{i}=3150$ subjects who
carry the variant rs12916-T. In this dataset, the sample size in the
controlled group $n_{0}$ is less than the dimension $p$. Our response variable
$Y_{i}$ is the individual plasma LDL-cholesterol concentration measured in
mg/dl.

As for the covariates $X_{i}$ for the subject $i$, we include the following
variables: $X_{i1}$ is the intercept, $X_{i2}$ represents age, $X_{i3}$
represents number of days of moderate physical activity, $X_{i4}$ represents
number of days of vigorous physical activity, $X_{i5}$ represents cooked
vegetable intake, $X_{i6}$ represents salad/raw vegetable intake, $X_{i7}$
represents fresh fruit intake, $X_{i8}$ represents dried fruit intake,
$X_{i9}$ represents oily fish intake, $X_{i,10}$ represents non-oily fish
intake, $X_{i,11}$ represents processed meat intake, $X_{i,12}$ represents
poultry intake, $X_{i,13}$ represents beef intake, $X_{i,14}$ represents
lamb/mutton intake, $X_{i,15}$ represents pork intake, $X_{i,16}$ represents
alcohol intake frequency per week, $X_{i,17}$ represents smoking status,
$X_{i,19}$ represents insulin medication usage, $X_{i,18}$ represents gender,
and $X_{i,20},\ldots,X_{i,p}$ contain additional 619 SNPs associated with the
LDL cholesterol concentration. The unit measurement of the included dietary
variables is tablespoons/day. We have provided detailed data pre-processing
steps in the Supplementary Materials.

Since our goal is to investigate whether statin usage has differential effects
on the study cohort, we estimate the HQTE
$\alpha(\tau;z)=z^{\prime}(\theta_{1}(\tau)-\theta_{0}(\tau))$ for two
different sets of the vector $z$:

In the first design, we study whether the combined effect of healthy dietary
patterns, physical activities, and reduced smoking differs in the statin-
treated and control groups for lowering LDL levels. We thus set

$\displaystyle
z=(0,0,\underbrace{1,\ldots,1}_{8},\underbrace{-1,\ldots,-1}_{6},0\ldots,0)^{\prime}\in\mathbb{R}^{p}.$

Figure 3(A) shows the estimated linear combination of quantile regression
coefficients $z^{\prime}\hat{\theta}_{1}(\tau)$ (green curve) and
$z^{\prime}\hat{\theta}_{0}(\tau)$ (red curve) along with estimated uniform
95% confidence bands (see Example 7 for details on how they were constructed).
We observe that the effect of statin usage, moderate to vigorous physical
activity combined with healthier dietary patterns on reducing the plasma LDL
cholesterol concentration is the largest among those patients whose
cholesterol are in the upper quantiles (i.e. whose cholesterol levels are high
relative to the population). For subjects who do not take statins, the effect
of increased moderate and vigorous physical activity combined with healthier
dietary patterns on reducing the plasma LDL cholesterol concentration is
roughly a quadratic function of the quantile $\tau$, but the effect overall
seems to be marginal. Figure 3(B) shows the estimate of $\alpha(\tau;z)$. We
see that the effect of statin usage is heterogeneous across different
quantiles of the LDL cholesterol concentration – its influence is more
significant at the right tail of the distribution. This suggests that statin
usage may help further reduce the LDL cholesterol level when combined with
healthy lifestyles for AD patients with rather high LDL cholesterol
concentration. Our findings might be helpful for researchers to design future
clinical trials to study the effect of statin usage on patients with high LDL
concentrations at baseline.

Figure 4: Heterogeneous quantile treatment effects of statin usage for three
subjects in the considered study sample with AD from the UK Biobank study.
Uniform 95% confidence bands discussed in Example 7 are given by the shaded
regions.

In the second design, we estimate HQTEs for three study participants in the UK
Biobank cohorts with different lifestyles: The first patient is a 65-year-old
subject (A) who exercises 4 times a week, has sufficient cooked and raw
vegetable intake (4 tablespoons per day), has one tablespoon of fruit intake
per day, has 2 drinks per week and with a recent smoking history. The second
patient is a 70-year-old subject (B) who exercises every day, has sufficient
vegetable and fruit intake (more than 3 tablespoons per day), has 2 drinks per
week and with a recent smoking history. The third patient is a 66-year-old
subject (C) who has no moderate/vigorous physical activity, less than 2
tablespoons vegetable and fruit intake per day, 5 drinks per week and no
recent smoking history. We do not observe a notable difference for meat intake
among these three subjects. The point estimates and their corresponding 95%
uniform confidence bands are reported in Figure 4. We have excluded these
three individuals from implementing our rank-score debiasing procedure.
Although the observed differential effects across the above three subjects
might be ascribed to the fact that study participants have different
genotypes, our study results suggest that the benefit of statin usage can be
heterogeneous across study participants.

##  8 Discussion

In this article, we have introduced a new procedure to study treatment effect
heterogeneity based on quantile regression modeling and rank-score debiasing.
While our rank-score debiased estimator is easy to implement and enjoys strong
theoretical guarantees, the following points merit future research:

First, it is worthwhile to relax the unconfoundedness assumption, simply
because unmeasured confounding presents a critical challenge to causal
inference from observational studies. A classical approach to mitigate the
confounding bias is to include instrumental variable methods. In this context,
the identification condition 1 can be modified similar to that of Chernozhukov
and Hansen, (2005). In future work, we therefore plan to investigate how to
combine instrumental variables with our proposed debiasing procedure.

Second, it is desirable to study the asymptotic efficiency of the rank-score
debiasing procedure. Since the rank-score debiased estimate is a scalar-valued
functional of a high-dimensional (nuisance) parameter, this requires a semi-
parametric point of view. However, the existing semiparametric efficiency
bounds for quantile regression apply only to fixed dimensional settings when
the quantile regression vector is independent of the sample size (Newey and
Powell,, 1990; Zhao,, 2001). The treatment of high-dimensional quantile
regression models requires a more elaborate analysis since the quantile
regression vector may change with sample size. In future work we intend to
develop a useful concept of semi-parametric efficiency of high-dimensional
processes indexed by changing function classes to answer the question whether
the rank-score debiased estimator is also asymptotically optimal in high
dimensions.

### Acknowledgement

We are deeply grateful for the associate editor and two anonymous reviewers
for their constructive suggestions that helped to improve presentation and
technical accuracy of this work. We thank the individuals involved in the UK
Biobank for their participation and the research teams for their work on
collecting, processing, and sharing these datasets. This research has been
conducted using the UK Biobank Resource (application number 48240), subject to
a data transfer agreement. Alexander Giessing thanks Jianqing Fan for his
unconditional support as Postdoctoral Research Associate at the Department of
Operations Research and Financial Engineering at Princeton University.
Jingshen Wang also acknowledges the support of the National Science Foundation
(DMS 2015325) and the National Institute of Health (R01MH125746).

## References

  * Abadie et al., (2002) Abadie, A., Angrist, J., and Imbens, G. (2002).  Instrumental variables estimates of the effect of subsidized training on the quantiles of trainee earnings.  Econometrica, 70(1):91–117.
  * Adamczak, (2008) Adamczak, R. (2008).  A tail inequality for suprema of unbounded empirical processes with applications to markov chains.  Electronic Journal of Probability, 13:1000–1034.
  * Alzheimer’s Association, (2019) Alzheimer’s Association, A. A. (2019).  2019 alzheimer’s disease facts and figures.  Alzheimer’s & dementia, 15(3):321–387.
  * Angrist et al., (2006) Angrist, J., Chernozhukov, V., and Fernández-Val, I. (2006).  Quantile regression under misspecification, with an application to the U.S. wage structure.  Econometrica, 74(2):539–563.
  * Angrist, (2004) Angrist, J. D. (2004).  Treatment effect heterogeneity in theory and practice.  The economic journal, 114(494):C52–C83.
  * Athey et al., (2018) Athey, S., Imbens, G. W., and Wager, S. (2018).  Approximate residual balancing: debiased inference of average treatment effects in high dimensions.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(4):597–623.
  * Atkins et al., (2017) Atkins, Summer, Einarsson, Gudmundur, Ames, Brendan, Clemmensen, and Line (2017).  Proximal methods for sparse optimal scoring and discriminant analysis.  arXiv preprint arXiv:1705.07194.
  * Bahadur, (1966) Bahadur, R. R. (1966).  A Note on Quantiles in Large Samples.  The Annals of Mathematical Statistics, 37(3):577 – 580.
  * Barthold et al., (2020) Barthold, D., Joyce, G., Diaz Brinton, R., Wharton, W., Kehoe, P. G., and Zissimopoulos, J. (2020).  Association of combination statin and antihypertensive therapy with reduced alzheimer’s disease and related dementia risk.  PloS one, 15(3):e0229541.
  * Belloni and Chernozhukov, (2011) Belloni, A. and Chernozhukov, V. (2011).  $l_{1}$-penalized quantile regression in high-dimensional sparse models.  The Annals of Statistics, 39(1):82–130.
  * Belloni and Chernozhukov, (2013) Belloni, A. and Chernozhukov, V. (2013).  Least squares after model selection in high-dimensional sparse models.  Bernoulli, 19(2):521–547.
  * (12) Belloni, A., Chernozhukov, V., Chetverikov, D., and Fernández-Val, I. (2019a).  Conditional quantile processes based on series or many regressors.  Journal of Econometrics, 213(1):4 – 29.  Annals: In Honor of Roger Koenker.
  * (13) Belloni, A., Chernozhukov, V., and Kato, K. (2019b).  Valid post-selection inference in high-dimensional approximately sparse quantile regression models.  Journal of the American Statistical Association, 114(526):749–758.
  * Bickel et al., (2009) Bickel, P. J., Ritov, Y., and Tsybakov, A. B. (2009).  Simultaneous analysis of lasso and dantzig selector.  The Annals of Statistics, 37(4):1705–1732.
  * Bowers et al., (2011) Bowers, K., Liu, G., Wang, P., Ye, T., Tian, Z., Liu, E., Yu, Z., Yang, X., Klebanoff, M., and Yeung, E. (2011).  Birth weight, postnatal weight change, and risk for high blood pressure among chinese children.  Pediatrics, 127(5):e1272–e1279.
  * Bradic and Kolar, (2017) Bradic, J. and Kolar, M. (2017).  Uniform inference for high-dimensional quantile regression: linear functionals and regression rank scores.  arXiv preprint arXiv:1702.06209.
  * Breiman et al., (1984) Breiman, L., Friedman, J., Stone, C. J., and Olshen, R. A. (1984).  Classification and regression trees.  CRC press.
  * Cattaneo, (2010) Cattaneo, M. D. (2010).  Efficient semiparametric estimation of multi-valued treatment effects under ignorability.  Journal of Econometrics, 155(2):138–154.
  * Chao et al., (2017) Chao, S.-K., Volgushev, S., and Cheng, G. (2017).  Quantile processes for semi and nonparametric regression.  Electronic Journal of Statistics, 11(2):3272–3331.
  * Chen, (2007) Chen, X. (2007).  Chapter 76 large sample sieve estimation of semi-nonparametric models.  volume 6 of Handbook of Econometrics, pages 5549–5632. Elsevier.
  * Chernozhukov et al., (2018) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2018).  Double/debiased machine learning for treatment and structural parameters: Double/debiased machine learning.  The Econometrics Journal, 21(1).
  * Chernozhukov et al., (2014) Chernozhukov, V., Chetverikov, D., and Kato, K. (2014).  Gaussian approximation of suprema of empirical processes.  The Annals of Statistics, 42(4):1564–1597.
  * Chernozhukov and Fernández-Val, (2005) Chernozhukov, V. and Fernández-Val, I. (2005).  Subsampling inference on quantile regression processes.  Sankhya: The Indian Journal of Statistics (2003-2007), 67(2):253–276.
  * Chernozhukov and Hansen, (2005) Chernozhukov, V. and Hansen, C. (2005).  An iv model of quantile treatment effects.  Econometrica, 73(1):245–261.
  * Coppock et al., (2018) Coppock, A., Leeper, T. J., and Mullinix, K. J. (2018).  Generalizability of heterogeneous treatment effect estimates across samples.  Proceedings of the National Academy of Sciences, 115(49):12441–12446.
  * Di Paolo and Kim, (2011) Di Paolo, G. and Kim, T.-W. (2011).  Linking lipids to alzheimer’s disease: cholesterol and beyond.  Nature Reviews Neuroscience, 12(5):284–296.
  * Dudley, (2014) Dudley, R. (2014).  Uniform Central Limit Theorems.  Cambridge Studies in Advanced Mathematics. Cambridge University Press.
  * Fan and Li, (2001) Fan, J. and Li, R. (2001).  Variable selection via nonconcave penalized likelihood and its oracle properties.  Journal of the American statistical Association, 96(456):1348–1360.
  * Firpo, (2007) Firpo, S. (2007).  Efficient semiparametric estimation of quantile treatment effects.  Econometrica, 75(1):259–276.
  * Foucart and Rauhut, (2013) Foucart, S. and Rauhut, H. (2013).  A Mathematical Introduction to Compressive Sensing.  Springer New York.
  * Frölich and Melly, (2013) Frölich, M. and Melly, B. (2013).  Unconditional quantile treatment effects under endogeneity.  Journal of Business & Economic Statistics, 31(3):346–357.
  * Fu et al., (2017) Fu, A., Narasimhan, B., and Boyd, S. (2017).  Cvxr: An r package for disciplined convex optimization.  arXiv preprint arXiv:1711.07582.
  * Geifman et al., (2017) Geifman, N., Brinton, R. D., Kennedy, R. E., Schneider, L. S., and Butte, A. J. (2017).  Evidence for benefit of statins to modify cognitive decline and risk in alzheimer’s disease.  Alzheimer’s research & therapy, 9(1):1–10.
  * Giessing, (2022) Giessing, A. (2022).  Inequalities for suprema of unbounded empirical processes.  Technical Report.
  * Giné and Nickl, (2015) Giné, E. and Nickl, R. (2015).  Mathematical Foundations of Infinite-Dimensional Statistical Models.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press.
  * He and Shi, (1994) He, X. and Shi, P. (1994).  Convergence rate of b-spline estimators of nonparametric conditional quantile functions.  Journal of Nonparametric Statistics, 3(3-4):299–308.
  * He et al., (2013) He, X., Wang, L., and Hong, H. G. (2013).  Quantile-adaptive model-free variable screening for high-dimensional heterogeneous data.  Annals of Statistics, 41(1):342–369.
  * Imai and Ratkovic, (2013) Imai, K. and Ratkovic, M. (2013).  Estimating treatment effect heterogeneity in randomized program evaluation.  The Annals of Applied Statistics, 7(1):443–470.
  * Jack Jr et al., (2010) Jack Jr, C. R., Knopman, D. S., Jagust, W. J., Shaw, L. M., Aisen, P. S., Weiner, M. W., Petersen, R. C., and Trojanowski, J. Q. (2010).  Hypothetical model of dynamic biomarkers of the alzheimer’s pathological cascade.  The Lancet Neurology, 9(1):119–128.
  * Jick et al., (2000) Jick, H., Zornberg, G. L., Jick, S. S., Seshadri, S., and Drachman, D. A. (2000).  Statins and the risk of dementia.  The Lancet, 356(9242):1627–1631.
  * Kern et al., (2016) Kern, H. L., Stuart, E. A., Hill, J., and Green, D. P. (2016).  Assessing methods for generalizing experimental impact estimates to target populations.  Journal of research on educational effectiveness, 9(1):103–127.
  * Kiefer, (1967) Kiefer, J. (1967).  On Bahadur’s Representation of Sample Quantiles.  The Annals of Mathematical Statistics, 38(5):1323 – 1342.
  * Kivipelto et al., (2018) Kivipelto, M., Mangialasche, F., and Ngandu, T. (2018).  Lifestyle interventions to prevent cognitive impairment, dementia and alzheimer disease.  Nature Reviews Neurology, 14(11):653–666.
  * Kley et al., (2016) Kley, T., Volgushev, S., Dette, H., and Hallin, M. (2016).  Quantile spectral processes: Asymptotic analysis and inference.  Bernoulli, 22(3):1770 – 1807.
  * Koenker, (2005) Koenker, R. (2005).  Quantile Regression.  Econometric Society Monographs. Cambridge University Press, Cambridge.
  * Koenker and Zhao, (1994) Koenker, R. and Zhao, Q. (1994).  L-estimatton for linear heteroscedastic models.  Journaltitle of Nonparametric Statistics, 3(3-4):223–235.
  * Koltchinskii, (2011) Koltchinskii, V. (2011).  Oracle Inequalities in Empirical Risk Minimization and Sparse Recovery Problems: École D’Été de Probabilités de Saint-Flour XXXVIII-2008.  Lecture Notes in Mathematics. Springer Verlag, New York.
  * Künzel et al., (2019) Künzel, S. R., Sekhon, J. S., Bickel, P. J., and Yu, B. (2019).  Metalearners for estimating heterogeneous treatment effects using machine learning.  Proceedings of the national academy of sciences, 116(10):4156–4165.
  * Ledoux and Talagrand, (1996) Ledoux, M. and Talagrand, M. (1996).  Probability in Banach Spaces: Isoperimetry and Processes.  Springer-Verlag, Berlin.
  * Lipkovich et al., (2017) Lipkovich, I., Dmitrienko, A., and B D’Agostino Sr, R. (2017).  Tutorial in biostatistics: data-driven subgroup identification and analysis in clinical trials.  Statistics in medicine, 36(1):136–196.
  * Lourida et al., (2019) Lourida, I., Hannon, E., Littlejohns, T. J., Langa, K. M., Hyppönen, E., Kuźma, E., and Llewellyn, D. J. (2019).  Association of lifestyle and genetic risk with incidence of dementia.  Jama, 322(5):430–437.
  * Ma and Huang, (2017) Ma, S. and Huang, J. (2017).  A concave pairwise fusion approach to subgroup analysis.  Journal of the American Statistical Association, 112(517):410–423.
  * Maurer, (2016) Maurer, A. (2016).  A Vector-Contraction Inequality for Rademacher Complexities.  In Ortner, R., Simon, H. U., and Zilles, S., editors, Algorithmic Learning Theory, pages 3–17, Cham. Springer.
  * McGuinness et al., (2010) McGuinness, B., O’Hare, J., Craig, D., Bullock, R., Malouf, R., and Passmore, P. (2010).  Statins for the treatment of dementia.  Cochrane database of systematic reviews, (8).
  * Mhanna et al., (2015) Mhanna, M., Iqbal, A., and Kaelber, D. (2015).  Weight gain and hypertension at three years of age and older in extremely low birth weight infants.  Journal of neonatal-perinatal medicine, 8(4):363–369.
  * Newey and Powell, (1990) Newey, W. K. and Powell, J. L. (1990).  Efficient estimation of linear and type I censored regression models under conditional quantile restrictions.  Econometric Theory, pages 295–317.
  * Neyman, (1959) Neyman, J. (1959).  Optimal asymptotic tests of composite hypotheses.  Probability and statistics, pages 213–234.
  * Nie and Wager, (2019) Nie, X. and Wager, S. (2019).  Quasi-oracle estimation of heterogeneous treatment effects.  arXiv preprint arXiv:1712.04912.
  * Nissen et al., (2005) Nissen, S. E., Tuzcu, E. M., Schoenhagen, P., Crowe, T., Sasiela, W. J., Tsai, J., Orazem, J., Magorien, R. D., O’Shaughnessy, C., and Ganz, P. (2005).  Statin therapy, ldl cholesterol, c-reactive protein, and coronary artery disease.  New England Journal of Medicine, 352(1):29–38.
  * Panchenko, (2003) Panchenko, D. (2003).  Symmetrization approach to concentration inequalities for empirical processes.  The Annals of Probability, 31(4):2068–2081.
  * Pinelis, (2020) Pinelis, I. (2020).  Exact lower and upper bounds on the incomplete gamma function.
  * Rockwood et al., (2002) Rockwood, K., Kirkland, S., Hogan, D. B., MacKnight, C., Merry, H., Verreault, R., Wolfson, C., and McDowell, I. (2002).  Use of lipid-lowering agents, indication bias, and the risk of dementia in community-dwelling elderly people.  Archives of neurology, 59(2):223–227.
  * Rubin, (1974) Rubin, D. B. (1974).  Estimating causal effects of treatments in randomized and nonrandomized studies.  Journal of educational Psychology, 66(5):688.
  * Rubin, (2009) Rubin, D. B. (2009).  Should observational studies be designed to allow lack of balance in covariate distributions across treatment groups?  Statistics in Medicine, 28(9):1420–1423.
  * Semenova and Chernozhukov, (2021) Semenova, V. and Chernozhukov, V. (2021).  Debiased machine learning of conditional average treatment effects and other causal functions.  The Econometrics Journal, 24(2):264–289.
  * Swerdlow et al., (2015) Swerdlow, D. I., Preiss, D., Kuchenbaecker, K. B., Holmes, M. V., Engmann, J. E., Shah, T., Sofat, R., Stender, S., Johnson, P. C., Scott, R. A., et al. (2015).  Hmg-coenzyme a reductase inhibition, type 2 diabetes, and bodyweight: evidence from genetic analysis and randomised trials.  The Lancet, 385(9965):351–361.
  * van de Geer et al., (2014) van de Geer, S., Bühlmann, P., Ritov, Y., and Dezeure, R. (2014).  On asymptotically optimal confidence regions and tests for high-dimensional models.  The Annals of Statistics, 42(3):1166–1202.
  * van der Vaart and Wellner, (1996) van der Vaart, A. W. and Wellner, J. (1996).  Weak Convergence and Empirical Processes: With Applications to Statistics.  Springer Series in Statistics. Springer-Verlag, New York.
  * Wahlberg et al., (2012) Wahlberg, B., Boyd, S., Annergren, M., and Wang, Y. (2012).  An admm algorithm for a class of total variation regularized estimation problems.  IFAC Proceedings Volumes, 45(16):83–88.
  * Wang and He, (2021) Wang, L. and He, X. (2021).  Analysis of global and local optima of regularized quantile regression in high dimension: A subgradient approach.  Preprint.
  * Wang et al., (2018) Wang, L., Zhou, Y., Song, R., and Sherwood, B. (2018).  Quantile-optimal treatment regimes.  Journal of the American Statistical Association, 113(523):1243–1254.
  * Wang and Zubizarreta, (2017) Wang, Y. and Zubizarreta, J. (2017).  Minimal dispersion approximately balancing weights: Asymptotic properties and practical considerations.  Biometrika, 103(1):1–29.
  * WHO, (2019) WHO (2019).  Risk reduction of cognitive decline and dementia: Who guidelines.
  * Würtz et al., (2016) Würtz, P., Wang, Q., Soininen, P., Kangas, A. J., Fatemifar, G., Tynkkynen, T., Tiainen, M., Perola, M., Tillin, T., Hughes, A. D., et al. (2016).  Metabolomic profiling of statin use and genetic inhibition of hmg-coa reductase.  Journal of the American College of Cardiology, 67(10):1200–1210.
  * Zhang and Zhang, (2014) Zhang, C.-H. and Zhang, S. S. (2014).  Confidence intervals for low dimensional parameters in high dimensional linear models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76(1):217–242.
  * Zhao, (2001) Zhao, Q. (2001).  Asymptotically efficient median regression in the presence of heteroskedasticity of unknown form.  Econometric Theory, 17(4):765–784.
  * Zhao et al., (2019) Zhao, W., Zhang, F., and Lian, H. (2019).  Debiasing and distributed estimation for high-dimensional quantile regression.  IEEE Transactions on Neural Networks and Learning Systems.
  * Zhou and Portnoy, (1996) Zhou, K. Q. and Portnoy, S. L. (1996).  Direct use of regression quantiles to construct confidence sets in linear models.  Annals of Statistics, 24(1):287–306.
  * Zubizarreta, (2015) Zubizarreta, J. R. (2015).  Stable weights that balance covariates for estimation with incomplete outcome data.  Journal of the American Statistical Association, 110(511):910–922.

Supplementary Materials for “Debiased Inference on Heterogeneous Quantile
Treatment Effects with Regression Rank-Scores” Alexander
Giessing11footnotemark: 1 Jingshen Wang22footnotemark: 2

##  Appendix A Connection with existing literature

In this section, we briefly discuss the differences between our approach and
the ones propose by Belloni et al., 2019b  and Zhao et al., (2019). To
simplify the discussion, we introduce the following simple quantile regression
model: Let $\\{(Y_{i},X_{i})\\}_{i=1}^{n}$ be a random sample of response
variable $Y$ and covariates $X$. For a qantile level $\tau$, suppose that the
conditional quantile function of a response $Y$ conditional on the covariates
$X$ is given by

$\displaystyle
Q(\tau;x)=x^{\prime}\beta(\tau),\quad\quad{}\beta(\tau)=\big{(}\beta_{1}(\tau),\ldots,\beta_{p}(\tau)\big{)}^{\prime}\quad{}\quad{}\sup_{\tau\in\mathcal{T}}\left\|\beta(\tau)\right\|_{0}\ll
p\wedge n.$

We define the weighted covariance matrix

$\displaystyle
H(\tau)=\mathbb{E}\big{[}f_{i}^{2}(\tau)X_{i}X_{i}^{\prime}\big{]}=\mathbb{E}\begin{bmatrix}f_{i}^{2}(\tau)X_{i1}^{2}&f_{i}^{2}(\tau)X_{i1}X_{i,-1}^{\prime}\\\
f_{i}^{2}(\tau)X_{i1}X_{i,-1}&f_{i}^{2}(\tau)X_{i,-1}X_{i,-1}^{\prime}\end{bmatrix}\equiv\begin{bmatrix}H_{11}(\tau)&H_{12}(\tau)\\\
H_{21}(\tau)&H_{22}(\tau)\end{bmatrix},$

where $f_{i}(\tau)$ is the conditional density of $Y|X$ evaluated at
$X_{i}^{\prime}\beta(\tau)$.

Belloni et al., 2019b  provide a procedure for inference on a single
regression coefficient, say $\beta_{1}(\tau)$. When $p$ is larger than the
sample size $n$ and model selection is adopted to carry out inference on
$\beta_{1}(\tau)$, the authors propose a method that is robust to model
selection mistakes by constructing orthogonal score functions. They
hypothesize the auxiliary regression problem

$\displaystyle
f_{i}(\tau)X_{i1}=f_{i}(\tau)X_{i,-1}^{\prime}\gamma(\tau)+\nu_{i},\quad{}\quad{}\mathbb{E}[f_{i}(\tau)X_{i1}\nu_{i}]=0,\quad{}\quad{}i=1,\ldots,n,$
$\displaystyle\gamma(\tau)=\arg\min_{\gamma}\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}[f_{i}^{2}(\tau)(X_{i1}-X_{i,-1}^{\prime}\gamma)^{2}],$

where $X_{i,-1}=(X_{i2},\ldots,X_{ip})^{\prime}$ with regression vector
$\gamma(\tau)=H_{22}^{-1}(\tau)H_{21}(\tau)$ and residuals
$\nu_{i}=f_{i}(\tau)X_{i1}-f_{i}(\tau)X_{i,-1}^{\prime}H_{22}^{-1}(\tau)H_{21}(\tau)$.
This auxiliary regression motivates the authors to define the orthogonal score
function $\psi_{\tau}(Y,X^{\prime}b,\nu)=(\tau-\mathbf{1}\\{Y\leq
X^{\prime}b\\})\nu$. This score function not only satisfies
$\mathbb{E}[\psi_{\tau}(Y_{i},X_{i}^{\prime}\beta(\tau),\nu_{i})]=0$ but (by
the very definition of the regression vector $\gamma(\tau)$ and the residuals
$\nu_{i}$) also the orthogonality conditions

$\displaystyle\frac{\partial}{\partial
b}\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}\big{[}\big{(}\tau-\mathbf{1}\\{Y_{i}\leq
X_{i,1}\beta_{1}(\tau)+X_{i,-1}^{\prime}b\\}\big{)}\cdot\nu_{i}\big{]}\bigg{|}_{b=\beta_{-1}(\tau)}=0,$
$\displaystyle\frac{\partial}{\partial\gamma}\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}\big{[}(\tau-\mathbf{1}\\{Y_{i}\leq
X_{i}^{\prime}\beta(\tau)\\}\big{)}\cdot
f_{i}(\tau)\cdot\big{(}X_{i1}-X_{i,-1}^{\prime}\gamma\big{)}\big{]}\bigg{|}_{\gamma=\gamma(\tau)}=0.$

The authors proceed to use this score function to construct an estimator
$\hat{\beta}_{1}(\tau)$ which, by virtue of above orthogonality conditions, is
protected against the estimation biases of the high-dimensional nuisance
parameters $\gamma(\tau)$ and
$\big{(}\beta_{2}(\tau),\ldots,\beta_{p}(\tau)\big{)}^{\prime}$ (Belloni et
al., 2019b, , Section 2). The theoretical investigation of Belloni et al.,
2019b  shows that the asymptotic variance of their estimator equals

$\displaystyle\tau(1-\tau)\bar{\mathbb{E}}[\nu_{i}^{2}]^{-1}=\tau(1-\tau)\cdot\big{[}H_{11}(\tau)-H_{12}(\tau)H_{22}^{-1}(\tau)H_{21}(\tau)\big{]}^{-1}=\tau(1-\tau)[H(\tau)^{-1}]_{11}.$

Comparing this result with our Theorem 1, we see that when the goal is to
estimate a single regression coefficient, our rank-score debiased estimator is
asymptotically equivalent to the one proposed by Belloni et al., 2019b . The
key difference between Belloni et al., 2019b  and our approach is that we
directly debias the linear combination of the high-dimensional quantile
regression coefficient $z^{\prime}\beta(\tau)$. We do not consider a low
dimensional projection of the high-dimensional regression vector.

Zhao et al., (2019) propose a debiased penalized quantile regression
coefficient estimator by directly estimating the weighted inverse covariance
matrix $\mathbb{E}\big{[}f_{i}(\tau)X_{i}X_{i}^{\prime}\big{]}$ based on ideas
from van de Geer et al., (2014) and Zhang and Zhang, (2014). Concretely, their
estimator for $\beta(\tau)$ is of the form

$\displaystyle\hat{\beta}_{\mathrm{Lasso}}(\tau)=\arg\min_{b\in\mathbb{R}^{p}}\left\\{\frac{1}{n}\sum_{i=1}^{n}\rho_{\tau}\big{(}Y_{i}-X_{i}^{\prime}b\big{)}+p_{\lambda}(b)\right\\},$
$\displaystyle\hat{\beta}(\tau)=\hat{\beta}_{\mathrm{Lasso}}(\tau)+\widehat{\Theta}\left(\frac{1}{n}\sum_{i=1}^{n}\big{(}\tau-\mathbf{1}\\{Y_{i}\leq
X_{i}^{\prime}\hat{\beta}_{\texttt{Lasso}}(\tau)\\}\big{)}X_{i}\right),$

where $\widehat{\Theta}$ is an estimate of the inverse covariance matrix
$D_{1}(\tau)^{-1}:=\Big{[}\mathbb{E}\big{[}f_{i}(\tau)X_{i}X_{i}^{\prime}\big{]}\Big{]}^{-1}$
(the authors use the R package clime to estimate the inverse covariance
matrix). The authors establish consistency of this estimator in the
$\ell_{\infty}$-norm but do not derive weak convergence results. Based on the
theoretical derivation given in their Appendix A, their estimator satisfies
that

$\displaystyle\hat{\beta}_{1}(\tau)-\beta_{1}(\tau)=e_{1}^{\prime}\widehat{\Theta}\left(\frac{1}{n}\sum_{i=1}^{n}\big{(}\tau-\mathbf{1}\\{Y_{i}\leq
X_{i}^{\prime}\hat{\beta}_{\texttt{Lasso}}(\tau)\\}\big{)}X_{i}\right)+r_{n},$

where $e_{1}=(1,0,\ldots,0)^{\prime}$ and $r_{n}$ is a remainder term such
that $\|r_{n}\|_{\infty}\downarrow 0$. The rate at which $\|r_{n}\|_{\infty}$
vanishes is too slow for any weak convergence guarantees, but if it was fast
enough, the asymptotic variance of their estimator would be of the typical
sandwich-formula

$\displaystyle\tau(1-\tau)\cdot
e_{1}^{\prime}D_{1}^{-1}(\tau)D_{0}D_{1}^{-1}(\tau)e_{1},$

where $D_{0}:=\mathbb{E}\big{[}X_{i}X_{i}^{\prime}\big{]}$. This indicates
that the estimator in Zhao et al., (2019) is not as efficient as the rank-
score debiased estimator (see also Koenker,, 2005, Section 5.3.1). As an
aside, this asymptotic variance can obviously be estimated as
$\frac{1}{n}\sum_{i=1}^{n}(e_{1}^{\prime}\widehat{\Theta}X_{i})^{2}$.

##  Appendix B Additional implementation details for the case study

We first provide a pipeline for cleaning the UK Biobank data for our case
study. We obtain phenotype and genotype data from the UK Biobank study with
following steps. For the phenotype data, we download encoded data in .enc
format from the UK Biobank’s Access Management System (AMS). To decrypt the
encoded data, we download three helper programs: ukb_md5, ukb_unpack, and
ukb_conv. Note that the helper programs are only supported by Windows and
Linux systems. Second, we verify the integrity of the encoded data via ukb_md5
and unpack them into .enc_ukb format with ukb_unpack. To convert the data into
a readible format, we use ukb_conv to convert the .enc_ukb data into .csv
format (other options include txt, r, sas stata or bulk format). The data
dictionary can be obtained using ukb_conv with the docs option. After
decrypting and converting the encoded data, we obtain a dataset with sample
size n = $502,481$ and $20,502$ features. In our study, we work with the
phenotypes at the baseline to avoid confounding issues, and we extract
baseline variables with the suffix “-0.0”. The phenotype data we extract
include gender, age at recruitment, AD family history, International
Statistical Classification of Disease 9th revision (ICD 9) and 10th revision
(ICD 10) codes and self-report for T2D and AD.

For the genotype data, we download imputed genotypes and associated sample
information for 23 chromosomes. Imputation BGEN and Imputation sample can be
obtained via ukb_gene program. Imputation BGI and Imputation MAF+info can be
downloaded directly from UK Biobank resources $1965$ and $1967$. Second, we
use the Imputation sample file to remove individuals without genotype
information, which yields sample size $n=407,057$. Finally, we read in BGEN
files with snp_readBGEN function in R package bigsnpr using HapMap3 as the
reference genomes. snp_readBGEN converts the BGEN files into an R object
comprising of two elements: genotype and map, where genotype represents the
imputed genotypes in a matrix format and map contains the features of SNPs
(chromosome, rsid, physical position, major and minor alleles and allele
frequency). We only extract the genotype matrix from genotype and rsid from
map as our genotype data. We restrict our sample to subjects used for the
principle components (PCs) computation, since those individuals are unrelated.
From the extracted genotype data, we obtain the treatment variable: rs12916
(on chromosome 5), a functionally equivalent SNP of statins.

Since rs12916 and LDL cholestrol concertration are associated with other SNPs
due to genetic pleiotropy and linkage disequilibrium, we adjust for low-
density lipoprotein (LDL) and rs12916 related SNPs in our study. To find
disease-associated SNPs, we rely on published GWAS studies from the GWAS
catalogue. In our study, we define the disease-associated SNPs as SNPs
associated with LDL with $p$-values less than $5\times 10^{-5}$. To determine
the $p-$values for multiple correlated SNPs in the same locus, we use the
linkage disequilibrium clumping procedure with $R^{2}<0.01$. Our filtration
criteria yield 619 disease-associated SNPs.

Lastly, we note that as the R package quantreg with the “Lasso” option tends
to be numerically unstable in our real data analyses (it produces error
message “singular design matrix”), we switch to the the R package conquer for
the implementation of the quantile regression problems.

##  Appendix C Overview of technical results

In the the remainder of these supplementary materials we develop the
theoretical backbone of the results in the main text. The key to the weak
convergence results in the main text is the Bahadur-type representation of the
rank-score balanced estimate of the CQF. Establishing this Bahadur-type
representation with a reasonable non-asymptotic bound on the remainder term is
very involved. In particular, the proof requires the following auxiliary
results:

  * •

uniform consistency of the _unweighted_ $\ell_{1}$-penalized quantile
regression vector;

  * •

uniform control over the empirical sparsity of the _unweighted_
$\ell_{1}$-penalized quantile regression vector;

  * •

uniform consistency of the solution to the dual program of the rank-score
balancing program.

These results are new. Existing results and proofs concerning uniform
consistency and empirical sparsity of the _weighted_ $\ell_{1}$-penalized
quantile regression vector (Belloni and Chernozhukov,, 2011; Belloni et al.,