direction switch/.style= for tree=edge+=thick, font=, where level>=1folder,
grow’=0for children=forked edge, where level=3draw, ,

# Assumption-Lean Quantile Regression

Georgi Baklicharov
Department of Applied Mathematics, Computer Science and Statistics
Ghent University
Krijgslaan 281 (S9), 9000 Ghent, Belgium
<EMAIL_ADDRESS>&Christophe Ley
Department of Mathematics
Université du Luxembourg
Maison du Nombre, Avenue de la Fonte 6, L-4364 Esch-sur-Alzette, Luxembourg
<EMAIL_ADDRESS>&Vanessa Gorasso
Department of Epidemiology and Public Health
Sciensano
Rue J Wytsman 14, 1050 Brussels, Belgium
<EMAIL_ADDRESS>&Brecht Devleesschauwer
Department of Epidemiology and Public Health
Sciensano
Rue J Wytsman 14, 1050 Brussels, Belgium
<EMAIL_ADDRESS>&Stijn Vansteelandt
Department of Applied Mathematics, Computer Science and Statistics
Ghent University
Krijgslaan 281 (S9), 9000 Ghent, Belgium
<EMAIL_ADDRESS>
(April, 2024)

# Assumption-Lean Quantile Regression

Georgi Baklicharov
Department of Applied Mathematics, Computer Science and Statistics
Ghent University
Krijgslaan 281 (S9), 9000 Ghent, Belgium
<EMAIL_ADDRESS>&Christophe Ley
Department of Mathematics
Université du Luxembourg
Maison du Nombre, Avenue de la Fonte 6, L-4364 Esch-sur-Alzette, Luxembourg
<EMAIL_ADDRESS>&Vanessa Gorasso
Department of Epidemiology and Public Health
Sciensano
Rue J Wytsman 14, 1050 Brussels, Belgium
<EMAIL_ADDRESS>&Brecht Devleesschauwer
Department of Epidemiology and Public Health
Sciensano
Rue J Wytsman 14, 1050 Brussels, Belgium
<EMAIL_ADDRESS>&Stijn Vansteelandt
Department of Applied Mathematics, Computer Science and Statistics
Ghent University
Krijgslaan 281 (S9), 9000 Ghent, Belgium
<EMAIL_ADDRESS>
(April, 2024)

###### Abstract

Quantile regression is a powerful tool for detecting exposure-outcome
associations given covariates across different parts of the outcome’s
distribution, but has two major limitations when the aim is to infer the
effect of an exposure. Firstly, the exposure coefficient estimator may not
converge to a meaningful quantity when the model is misspecified, and
secondly, variable selection methods may induce bias and excess uncertainty,
rendering inferences biased and overly optimistic. In this paper, we address
these issues via partially linear quantile regression models which parametrize
the conditional association of interest, but do not restrict the association
with other covariates in the model. We propose consistent estimators for the
unknown model parameter by mapping it onto a nonparametric main effect
estimand that captures the (conditional) association of interest even when the
quantile model is misspecified. This estimand is estimated using the efficient
influence function under the nonparametric model, allowing for the
incorporation of data-adaptive procedures such as variable selection and
machine learning. Our approach provides a flexible and reliable method for
detecting associations that is robust to model misspecification and excess
uncertainty induced by variable selection methods. The proposal is illustrated
using simulation studies and data on annual health care costs associated with
excess body weight.

##  1 Introduction

Quantile regression, originally introduced by Koenker and Bassett (1978), is a
statistical method that provides a flexible framework for modeling the
relationship between a response variable and predictors. Unlike traditional
linear regression models that focus on estimating the conditional mean of the
response variable, quantile regression allows for the estimation of
conditional quantiles, which can provide valuable and broader insights into
the distributional characteristics of the response variable. It has gained
popularity in recent years thanks to its ability to capture heterogeneity in
the data and provide more robust estimates in the presence of outliers, skewed
and/or heavy-tailed distributions (Koenker, 2005, 2017; Koenker et al., 2017).
Quantile regression has been applied in a variety of fields, including
economics, finance, health, and environmental sciences, among others (e.g.
Koenker and Hallock, 2001; Yu et al., 2003) but deserves to be used more
frequently, because of the rich description it provides of various relevant
aspects of data distributions, combined with appealing ease of interpretation
(Das et al., 2019).

In this work, we focus on the estimation of causal quantile treatment or
exposure effects under standard causal assumptions (Hernán and Robins, 2020).
The quantile treatment effect has originally been defined by Doksum (1974) and
Lehmann and D’Abrera (1975) in a two-sample setting as the horizontal distance
(i.e. the distance along the horizontal axis between two fixed points of
different distribution functions) between the outcome distribution of exposed
versus unexposed. However, extension to allow for covariate adjustment via
standard quantile regression raises concerns due to potential model
misspecification and post-selection inference.

Model misspecification is a concern for most statistical methods. When the
goal is solely prediction, this is often less problematic. Wrong models can
still deliver good predictions, because predictions are commonly obtained by
minimizing a measure of prediction error over the considered model class.
Moreover, unbiased estimates of prediction error are obtained, even when the
model is misspecified. However, when estimating exposure effects, model
misspecification typically causes more severe problems. For instance, in
(generalized) linear regression models that have a main exposure effect and
covariates, the ordinary least squares estimator generally fails to converge
to a comprehensive measure of the conditional relationship between the
exposure and outcome in the presence of model misspecification (Vansteelandt
and Dukes, 2022). In quantile regression models, the situation is even worse
as there is no good understanding of where the effect of a randomized
treatment converges to when the model is wrong. One possible solution is to
make the quantile model sufficiently complex by including interactions,
higher-order terms, transformations, etc. However, it is generally unclear how
much additional complexity is needed, and moreover, this added complexity may
lead to overfitting and thus estimates with large bias and standard errors. In
contrast, simple models help to make results insightful and easy to
communicate. In practice, one is therefore often torn between keeping the
model simple so that it is comprehensible and making it complex so that it
represents the ground truth (Breiman, 2001). A common way of finding a trade-
off between model simplicity and model correctness is by using data-adaptive
procedures, such as variable selection algorithms. However, this also raises
concerns because data-adaptive procedures usually give rise to additional
uncertainty that is not taken into account in standard error calculations. In
particular, standard inference commonly ignores the fact that making decisions
on whether or not to add certain covariates was data-driven, while it is
highly unlikely to make the exact same decisions under repeated sampling (Leeb
and Pötscher, 2005). This problem of post-selection inference renders standard
inferences overly optimistic.

In previous literature, attempts have been made to circumvent this issue by
employing nonparametric approaches, focusing on marginal effects. However, the
majority of these approaches is limited to binary exposures (e.g. Melly, 2006;
Firpo, 2007; Chernozhukov et al., 2013; Donald and Hsu, 2014; Caracciolo and
Furno, 2017; Athey et al., 2023). Ai et al. (2022) extended the work of Donald
and Hsu (2014) to continuous exposures, but still focussing on marginal
quantile effects. Instead, our objective is to infer conditional quantile
effects as these are less vulnerable to extrapolation (see the discussion
section). Our proposal will be in line with recent developments in causal
machine learning (Díaz, 2017; Belloni et al., 2017; Kallus et al., 2019) in
the sense that we allow the use of machine learning techniques, while ensuring
valid post-selection inference.

Throughout this work, we will focus on the conditional association between an
outcome $Y$ and a possibly continuous exposure $A$, given a vector of measured
covariates denoted by $L$ that are sufficient to control for confounding.
Therefore, we consider the partially linear quantile model

$\displaystyle Q_{\tau}(Y\mid A,L)=\beta_{\tau}A+\omega_{\tau}(L),$ (1)

where $Q_{\tau}(Y\mid A,L)$ denotes the $\tau$-quantile of the conditional
distribution of $Y$ given $A$ and $L$ for a given value of $0<\tau<1$,
$\beta_{\tau}$ is an unknown scalar parameter and $\omega_{\tau}(.)$ is an
unknown function. Both $\beta_{\tau}$ and $\omega_{\tau}(.)$ may be quantile-
specific. The parameter $\beta_{\tau}$ in this model reflects the conditional
association of interest, for which we will develop a debiased machine learning
estimator (Chernozhukov et al., 2018). First, we will assume model (1) to
hold. Later, we explain how our estimator still captures the association of
interest if the model is misspecified.

Our work builds on the large literature on partially linear quantile models
(e.g. Lee, 2003; Sun, 2005; Wu et al., 2010; Wu and Yu, 2014; Lv et al., 2015;
Sherwood and Wang, 2016; Zhong and Wang, 2023). Traditionally, two popular
approaches are being considered for estimating the nonparametric component of
the partially linear quantile model: kernel estimators (e.g. Lee, 2003; Sun,
2005; Lv et al., 2015), and splines (e.g. He and Shi, 1996; Sherwood and Wang,
2016). Those two approaches are also common in (partially linear) single-index
quantile models (e.g. Wu et al., 2010; Wu and Yu, 2014). However, those prior
developments have a number of important shortcomings. First, in high-
dimensional applications, they force high computational demands and challenges
stemming from reliance on kernel weighting or splines. To accommodate this
concern, Zhong and Wang (2023) proposed neural networks to estimate the
nonparametric component. However, a further concern with this and the
aforementioned proposals is that the probability limit of the estimators has
no straightforward interpretation when the partially linear quantile model is
misspecified (e.g., due to exposure-covariate interactions). In that case,
existing proposals also fail to provide valid data-adaptive inference (not
even for the probability limit of the estimator because it is a function of
high-dimensional nuisance parameters with complex distributions). Our proposal
provides a solution for both shortcomings by allowing the use of flexible
estimators that are known to work well in high dimensions, and by ensuring an
inference that is correct even under model misspecification. We take
inspiration from related work for generalized linear models in a discussion
paper by Vansteelandt and Dukes (2022) and for Cox models (Vansteelandt et
al., 2022), which are both closely related to recent developments on Targeted
Maximum Likelihood Estimation (TMLE) (van der Laan and Rose, 2011) and
Debiased Machine Learning (DML) (Chernozhukov et al., 2018).

##  2 Data-adaptive inference for $\beta_{\tau}$

###  2.1 A plug-in estimator

To construct a debiased machine learning estimator for $\beta_{\tau}$ in model
(1), we use the fact that it can be expressed as

$\Psi_{\tau}=\frac{E[(A-A^{*})\\{Q_{\tau}(Y\mid A,L)-Q_{\tau}(Y\mid
A^{*},L)\\}]}{E\\{(A-A^{*})^{2}\\}},$ (2)

when model (1) holds (see Appendix A for a proof). In this model-free
estimand, $A$ and $A^{*}$ are two exposure values independently drawn from
stratum $L$. In Section 3 we justify the choice of this estimand by arguing
that it is an interpretable measure of exposure effect, even if model (1) is
misspecified. Intuitively, this can already be seen because it is a weighted
average of differences in quantiles of the outcome for individuals with
different exposure values, but the same value for $L$; it is thus guaranteed
to summarize the conditional association between $A$ and $Y$ given $L$, even
if the model is wrong. For this nonparametric estimand, we will develop
nonparametric inference.

The model-free estimand $\Psi_{\tau}$ can also be written as (see Appendix A)

$\frac{E\left(\\{A-E(A\mid L)\\}\left[Q_{\tau}(Y\mid A,L)-E\\{Q_{\tau}(Y\mid
A,L)\mid L\\}\right]\right)}{E[\\{A-E(A\mid L)\\}^{2}]}.$ (3)

A so-called plug-in estimator is now straightforwardly obtained by
substituting $E(A\mid L)$, $Q_{\tau}(Y\mid A,L)$ and $E\\{Q_{\tau}(Y\mid
A,L)\mid L\\}$ by consistent estimators $\hat{E}(A\mid L)$,
$\hat{Q}_{\tau}(Y\mid A,L)$ and $\hat{E}\\{\hat{Q}_{\tau}(Y\mid A,L)\mid
L\\}$, respectively, and by replacing population expectations by sample
averages. This results in

$\frac{\sum_{i=1}^{n}\left(\\{A_{i}-\hat{E}(A_{i}\mid
L_{i})\\}\left[\hat{Q}_{\tau}(Y_{i}\mid
A_{i},L_{i})-\hat{E}\\{\hat{Q}_{\tau}(Y_{i}\mid A_{i},L_{i})\mid
L_{i}\\}\right]\right)}{\sum_{i=1}^{n}\\{A_{i}-\hat{E}(A_{i}\mid
L_{i})\\}^{2}}.$ (4)

To estimate the nuisance parameters, we will use flexible data-adaptive
procedures in line with recent developments on debiased machine learning (DML)
(van der Laan and Rose, 2011; Chernozhukov et al., 2018). These procedures may
include variable selection techniques for quantile regression, as well as more
advanced machine learning algorithms such as random forests or neural
networks.

However, plug-in estimators like (4) typically give rise to two problems. On
the one hand, data-adaptive procedures are optimized towards minimal
prediction error, but prediction is not our aim. Instead, our aim is to obtain
well-behaved estimators for this estimand. For example, a poor selection of
variables may induce confounding or selection bias. On the other hand, data-
adaptive estimators $\hat{E}(A\mid L)$, $\hat{Q}_{\tau}(Y\mid A,L)$ and
$\hat{E}\\{\hat{Q}_{\tau}(Y\mid A,L)\mid L\\}$ have a nonstandard
distribution, which is usually also unknown. This makes quantifying the
uncertainty in these estimators extremely difficult. Our proposal addresses
these concerns.

###  2.2 Debiased Machine Learning (DML)

In view of the aforementioned concerns, we aim to establish inference for our
main effect estimand (2) under the nonparametric model. This enables us to
avoid dependence on modeling constraints, which could result in excessively
optimistic inferences as well as computational difficulties. Like other
debiased machine learning strategies, we will make use of the efficient
influence function of the estimand to accomplish this goal (Pfanzagl, 1990;
Bickel et al., 1993). This is a mean zero functional of the observed data and
the data-generating distribution, which characterizes the estimand’s
sensitivity to arbitrary (smooth) changes in the data-generating law. Theorem
1, of which the proof can be found in Appendix B, gives the efficient
influence function of the estimand (3).

######  Theorem 1.

The efficient influence function of the main effect estimand $\Psi_{\tau}$,
defined by (3), under the nonparametric model is

$\displaystyle\begin{split}\phi(Y,A,L)&=\frac{A-E(A\mid L)}{E[\\{A-E(A\mid
L)\\}^{2}]}\bigg{[}Q_{\tau}(Y\mid A,L)-E\left\\{Q_{\tau}(Y\mid A,L)\mid
L\right\\}\\\ &\qquad\qquad+\frac{\tau-I\\{Y\leq Q_{\tau}(Y\mid
A,L)\\}}{f\\{Q_{\tau}(Y\mid A,L)\mid A,L\\}}-\Psi_{\tau}\\{A-E(A\mid
L)\\}\bigg{]},\end{split}$ (5)

where $I(.)$ and $f(.\mid A,L)$ denote an indicator function and the
conditional density of $Y$ given $A$ and $L$, respectively.

If the infinite-dimensional nuisance parameters $E(A\mid L)$, $Q_{\tau}(Y\mid
A,L)$, $E\left\\{Q_{\tau}(Y\mid A,L)\mid L\right\\}$ and $f\\{Q_{\tau}(Y\mid
A,L)\mid A,L\\}$ indexing the efficient influence function were known, then it
would follow from its mean zero property that a consistent estimator
$\tilde{\Psi}_{\tau}$ of $\Psi_{\tau}$ could be obtained as the value of
$\Psi_{\tau}$ that makes the sample average of the influence functions zero.
The resulting estimator’s asymptotic distribution would be governed by this
influence function in the sense that

$\displaystyle\sqrt{n}(\tilde{\Psi}_{\tau}-\Psi_{\tau})$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\frac{A_{i}-E(A_{i}\mid
L_{i})}{E[\\{A_{i}-E(A_{i}\mid L_{i})\\}^{2}]}\bigg{[}Q_{\tau}(Y_{i}\mid
A_{i},L_{i})-E\left\\{Q_{\tau}(Y_{i}\mid A_{i},L_{i})\mid L_{i}\right\\}$ (6)
$\displaystyle\qquad+\frac{\tau-I\\{Y_{i}\leq Q_{\tau}(Y_{i}\mid
A_{i},L_{i})\\}}{f\\{Q_{\tau}(Y_{i}\mid A_{i},L_{i})\mid
A_{i},L_{i}\\}}-\Psi_{\tau}\\{A_{i}-E(A_{i}\mid L_{i})\\}\bigg{]}+o_{p}(1).$
(7)

It follows that the estimator $\tilde{\Psi}_{\tau}$ is asymptotically normal
by the central limit theorem. Moreover, its bias approaches zero at a faster
rate than the standard error, and its variance can be estimated by computing
the sample variance of the efficient influence function. However, the
estimator $\tilde{\Psi}_{\tau}$ is practically infeasible because the
efficient influence function involves unknown nuisance parameters. We will
therefore substitute these with consistent data-adaptive estimators and denote
the resulting estimator $\hat{\Psi}_{\tau}$. In this way, we obtain the
closed-form estimator

$\displaystyle\begin{split}\hat{\Psi}_{\tau}&=\frac{1}{n}\sum_{i=1}^{n}\frac{A_{i}-\hat{E}(A_{i}\mid
L_{i})}{\frac{1}{n}\sum_{i=1}^{n}\\{A_{i}-\hat{E}(A_{i}\mid
L_{i})\\}^{2}}\Bigg{[}\hat{Q}_{\tau}(Y_{i}\mid A_{i},L_{i})\\\
&\qquad\qquad\qquad-\hat{E}\\{\hat{Q}_{\tau}(Y_{i}\mid A_{i},L_{i})\mid
L_{i}\\}+\frac{\tau-I\\{Y_{i}\leq\hat{Q}_{\tau}(Y_{i}\mid
A_{i},L_{i})\\}}{\hat{f}\\{\hat{Q}_{\tau}(Y_{i}\mid A_{i},L_{i})\mid
A_{i},L_{i}\\}}\Bigg{]}.\end{split}$ (8)

From the general theory on nonparametric estimation (Pfanzagl, 1990; Bickel et
al., 1993) and the results in the Appendix B, it follows that, when the
nuisance parameters are estimated on a sample independent to the one on which
$\hat{\Psi}_{\tau}$ is evaluated (to mitigate overfitting bias) and under
specific convergence rate assumptions detailed in Theorem 2, the estimator is
asymptotically linear in the sense that (6) holds with $\hat{\Psi}_{\tau}$
instead of $\tilde{\Psi}_{\tau}$. Consequently, its variance can be
consistently estimated as 1 over $n$ times the sample variance of the
influence functions, i.e. $1/n^{2}\sum_{i=1}^{n}\hat{\phi}(O_{i})^{2}$, where
$\hat{\phi}$ is defined like $\phi$ but with all nuisance parameters
substituted by consistent estimates.

Estimating the nuisance parameters on an independent sample can be attained
via sample splitting (Zheng and van der Laan, 2011) or cross-fitting
(Chernozhukov et al., 2018). In particular, the data can be split into $K$
disjoint subsets $I_{k},\,k=1,\ldots,K$. Let $I_{k}^{C}$ denote the complement
of the subset $I_{k}$. For each $k$, one can train data-adaptive nuisance
parameter estimators using only the data in $I_{k}^{C}$. In calculating
$\hat{\Psi}_{\tau}$, when a subject $i$ belongs to subset $I_{k}$, then the
nuisance parameter estimators trained on $I_{k}^{C}$ should be used to plug in
the nuisance parameter estimates in (8). When flexible data-adaptive estimates
are used for the nuisance parameters, then this cross-fitting is necessary for
the asymptotic validity of the proposal.

The advantage of utilizing the efficient influence function as the foundation
for the estimator is that it exhibits consistent asymptotic behavior,
regardless of whether it is based on the known nuisance parameters or
consistent estimates thereof. This is contingent on the convergence rate
conditions detailed in the following theorem. A proof is given in Appendix B.

######  Theorem 2.

The proposed estimator $\hat{\Psi}_{\tau}$ is asymptotically linear when

  1. 1.

all nuisance parameter estimators are consistent,

  2. 2.

all nuisance parameters are estimated on a sample independent to the one on
which $\hat{\Psi}_{\tau}$ is evaluated,

  3. 3.

all of the following terms are $o_{p}(n^{-1/2})$:

$\displaystyle E\left[\\{\hat{Q}_{\tau}(Y\mid A,L)-Q_{\tau}(Y\mid
A,L)\\}^{2}\right],$ (9) $\displaystyle
E\left(\left[1-\frac{f\\{Q_{\tau}(Y\mid A,L)\mid
A,L\\}}{\hat{f}\\{\hat{Q}_{\tau}(Y\mid A,L)\mid
A,L\\}}\right]^{2}\right)^{1/2}E\left[\\{\hat{Q}_{\tau}(Y\mid
A,L)-Q_{\tau}(Y\mid A,L)\\}^{2}\right]^{1/2},$ (10) $\displaystyle
E\left[\\{E(A\mid L)-\hat{E}(A\mid
L)\\}^{2}\right]^{1/2}E\left(\left[E(Q_{\tau}(Y\mid A,L)\mid
L)-\hat{E}\\{\hat{Q}_{\tau}(Y\mid A,L)\mid L\\}\right]^{2}\right)^{1/2},$
$\displaystyle E\left[\\{E(A\mid L)-\hat{E}(A\mid L)\\}^{2}\right],$

where the last condition is only needed when $\Psi_{\tau}\neq 0$, and where
the expectations are conditional on the training sample and thus only over the
distribution of $Y$, $A$ and $L$,

  4. 4.

$A-\hat{E}(A\mid L)=O_{p}(1)$,

  5. 5.

and $\;E[\\{A-E(A\mid L)\\}^{2}]>\epsilon$, $\;\hat{E}[\\{A-\hat{E}(A\mid
L)\\}^{2}]>\epsilon\;$ and $\;\hat{f}\\{\hat{Q}_{\tau}(Y\mid A,L)\mid
A,L\\}>\epsilon$ with probability 1 for some $\epsilon>0$.

The rate conditions stated in Theorem 2 will be satisfied e.g. if all nuisance
parameters are consistently estimated at a rate faster than $n^{1/4}$, but are
slightly less restrictive because of the rate conditions on product terms. The
third condition in Theorem 2 states that (9) needs to be $o_{p}(n^{-1/2})$.
However, if a faster convergence rate is achieved, then a slower convergence
rate for $\hat{f}\\{\hat{Q}_{\tau}(Y\mid A,L)\mid A,L\\}$ may suffice, as long
as the product (10) converges sufficiently fast. Similarly, if $E(A\mid L)$
can be estimated at a rate much faster than $n^{1/4}$, $E\\{Q_{\tau}(Y\mid
A,L)\mid L\\}$ can be estimated at a slower rate. Under particular
smoothness/sparsity assumptions, many data-adaptive methods can achieve these
convergence rates (a summary is given in Chernozhukov et al. (2018)).

The proposed estimator thus mitigates the plug-in bias that arises from the
bias-variance trade-off in flexible data-adaptive techniques. These techniques
are usually fine-tuned to minimize prediction error, leading to plug-in bias.
It is shown that (see e.g. van der Vaart (1998)) this plug-in bias can be
approximated as minus the sample average of the influence functions. Since our
estimator (8) is defined as the value of $\Psi_{\tau}$ that makes this sample
average zero, the plug-in bias will be approximately zero by construction. Our
proposal further takes into account the excess uncertainty induced by the use
of data-adaptive procedures (e.g. variable selection or machine learning
algorithms). By removing the bias induced by these procedures, our estimator
has the asymptotic behaviour of an oracle estimator where the nuisance
parameters are known. In this way, the proposed standard errors are honest in
that they reflect the uncertainty that arises from these data-driven model
choices.

###  2.3 Targeted Learning

Like other debiased machine learning strategies, targeted learning also
removes plug-in bias based on the efficient influence function, but it does so
after first ‘targeting’ the initial nuisance parameter estimates (van der Laan
and Rubin, 2006; van der Laan and Rose, 2011). It thereby delivers an
estimator that is asymptotically equivalent to the earlier defined debiased
estimator under the conditions of theorem 2. However, these estimators may
have different finite sample performance (Yadlowsky, 2022). Especially the
last term in (8) may sometimes be unstable because of the density in the
denominator. Therefore, we define a retargeted estimator
$\widetilde{Q}_{\tau}(Y\mid A,L)$ which forces

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left\\{A_{i}-\hat{E}(A_{i}\mid
L_{i})\right\\}\left[\frac{\tau-I\\{Y_{i}\leq\widetilde{Q}_{\tau}(Y_{i}\mid
A_{i},L_{i})\\}}{\hat{f}\\{\widetilde{Q}_{\tau}(Y_{i}\mid A_{i},L_{i})\mid
A_{i},L_{i}\\}}\right]=0.$ (11)

This is justified because (11) is an unbiased estimating equation for
$Q_{\tau}(Y\mid A,L)$. If (11) holds and $\hat{Q}_{\tau}$ is replaced by
$\widetilde{Q}_{\tau}$ in (8), the resulting estimator for $\Psi_{\tau}$
equals

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\frac{A_{i}-\hat{E}(A_{i}\mid
L_{i})}{\frac{1}{n}\sum_{i=1}^{n}\\{A_{i}-\hat{E}(A_{i}\mid
L_{i})\\}^{2}}\left[\widetilde{Q}_{\tau}(Y_{i}\mid
A_{i},L_{i})-\hat{E}\\{\widetilde{Q}_{\tau}(Y_{i}\mid A_{i},L_{i})\mid
L_{i}\\}\right],$

which is a plug-in estimator of estimand (3).

Constructing an estimator that satisfies (11) turns out to be rather
challenging. For convenience, we will write

$\displaystyle w(Q_{\tau,i})=\frac{A_{i}-\hat{E}(A_{i}\mid
L_{i})}{\hat{f}\\{Q_{\tau}(Y_{i}\mid A_{i},L_{i})\mid A_{i},L_{i}\\}}.$ (12)

Using the weights defined in (12), we can write equation (11) as

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}w(\widetilde{Q}_{\tau,i})\left[\tau-I\\{Y_{i}\leq\widetilde{Q}_{\tau}(Y_{i}\mid
A_{i},L_{i})\\}\right]=0.$

This equation needs to be solved for $\widetilde{Q}_{\tau}(Y_{i}\mid
A_{i},L_{i})$. Since it is highly nonlinear, we can only find an approximate
solution by updating the initial estimator $\hat{Q}_{\tau}(Y_{i}\mid
A_{i},L_{i})$ (e.g. obtained by using quantile regression forests) as follows:

  1. 1.

Let $Q_{\tau,i}^{(0)}=\hat{Q}_{\tau}(Y_{i}\mid A_{i},L_{i})$ and
$\epsilon_{0}=0$.

  2. 2.

For $k>0$: let
$Q_{\tau,i}^{(k)}=Q_{\tau,i}^{(k-1)}+\epsilon_{k-1}w\left(Q_{\tau,i}^{(k-1)}\right)$,
with

$\epsilon_{k-1}=\operatorname*{arg\,min}_{\epsilon}\left\lvert\frac{1}{n}\sum_{i=1}^{n}w\left(Q_{\tau,i}^{(k-1)}\right)\left[\tau-I\left\\{Y_{i}\leq
Q_{\tau,i}^{(k-1)}+\epsilon
w\left(Q_{\tau,i}^{(k-1)}\right)\right\\}\right]\right\rvert.$

  3. 3.

Repeat step 2 until the value of

$\left\lvert\frac{1}{n}\sum_{i=1}^{n}w\left(Q_{\tau,i}^{(k-1)}\right)\left[\tau-I\left\\{Y_{i}\leq
Q_{\tau,i}^{(k-1)}+\epsilon
w\left(Q_{\tau,i}^{(k-1)}\right)\right\\}\right]\right\rvert$

is no longer decreasing or achieves a certain threshold close to 0.

The resulting estimator for $Q_{\tau}(Y_{i}\mid A_{i},L_{i})$ is denoted
$\widetilde{Q}_{\tau}(Y_{i}\mid A_{i},L_{i})$. By now substituting
$\hat{Q}_{\tau}(Y_{i}\mid A_{i},L_{i})$ with $\widetilde{Q}_{\tau}(Y_{i}\mid
A_{i},L_{i})$ in the estimating equations estimator (8), we obtain a Targeted
Maximum Likelihood Estimation (TMLE) estimator

$\displaystyle\begin{split}\hat{\Psi}_{\tau}^{TMLE}&=\frac{1}{n}\sum_{i=1}^{n}\frac{A_{i}-\hat{E}(A_{i}\mid
L_{i})}{\frac{1}{n}\sum_{i=1}^{n}\\{A_{i}-\hat{E}(A_{i}\mid
L_{i})\\}^{2}}\left[\widetilde{Q}_{\tau}(Y_{i}\mid A_{i},L_{i})\right.\\\
&\left.\qquad\qquad-\hat{E}\\{\widetilde{Q}_{\tau}(Y_{i}\mid A_{i},L_{i})\mid
L_{i}\\}+\frac{\tau-I\\{Y_{i}\leq\widetilde{Q}_{\tau}(Y_{i}\mid
A_{i},L_{i})\\}}{\hat{f}\\{\widetilde{Q}_{\tau}(Y_{i}\mid A_{i},L_{i})\mid
A_{i},L_{i}\\}}\right].\end{split}$ (13)

The targeting step results in the last term of (13) being approximately zero,
such that

$\displaystyle\hat{\Psi}_{\tau}^{TMLE}\approx\frac{1}{n}\sum_{i=1}^{n}\frac{A_{i}-\hat{E}(A_{i}\mid
L_{i})}{\frac{1}{n}\sum_{i=1}^{n}\\{A_{i}-\hat{E}(A_{i}\mid
L_{i})\\}^{2}}\left[\widetilde{Q}_{\tau}(Y_{i}\mid
A_{i},L_{i})-\hat{E}\\{\widetilde{Q}_{\tau}(Y_{i}\mid A_{i},L_{i})\mid
L_{i}\\}\right].$ (14)

Comparing (14) with (4), it is clear that we removed (approximately) all bias
from the initial plug-in estimator by cleverly targeting the nuisance
parameter estimator for $Q_{\tau}(Y_{i}\mid A_{i},L_{i})$. In Appendix B, we
give a sketch of a proof that (13) and (14) are asymptotically equivalent,
even with only one iteration of targeting. Nevertheless, we opt to use (13) to
mitigate reliance on the accuracy of approximation (14).

The targeted learning estimator with 5-fold cross-fitting can be obtained as
follows. After splitting the data into 5 equal parts $I_{k}$, with
$k=1,\ldots,5$, five models are trained for each nuisance parameter, where
each time one subset $I_{k}$ is left out from the data. In particular, a
nonparametric model is trained for $E(A\mid L)$, $Q_{\tau}(Y\mid A,L)$,
$E\\{Q_{\tau}(Y\mid A,L)\mid L\\}$ and $f\\{Q_{\tau}(Y\mid A,L)\mid A,L\\}$ on
$I_{k}^{C}$ (i.e. 80% of the data). Then, those models are used to estimate
the nuisance parameters for subjects belonging to $I_{k}$. These estimates are
then used to calculate the contribution to $\Psi_{\tau}$ (which is a sample
average) for the same subjects. The targeting algorithm described above is
then applied to the entire data set. Then, $\Psi_{\tau}$ is estimated as in
(13). The standard error can be estimated similarly to the debiased estimator
(8), i.e. as the square root of

$\displaystyle\frac{1}{n^{2}}\sum_{i=1}^{n}\left(\frac{A_{i}-\hat{E}(A_{i}\mid
L_{i})}{\frac{1}{n}\sum_{i=1}^{n}\\{A_{i}-\hat{E}(A_{i}\mid
L_{i})\\}^{2}}\left[\widetilde{Q}_{\tau}(Y_{i}\mid
A_{i},L_{i})-\hat{E}\\{\widetilde{Q}_{\tau}(Y_{i}\mid A_{i},L_{i})\mid
L_{i}\\}\right.\right.$
$\displaystyle\left.\left.\qquad\qquad+\frac{\tau-I\\{Y_{i}\leq\widetilde{Q}_{\tau}(Y_{i}\mid
A_{i},L_{i})\\}}{\hat{f}\\{\widetilde{Q}_{\tau}(Y_{i}\mid A_{i},L_{i})\mid
A_{i},L_{i}\\}}-\hat{\Psi}_{\tau}^{TMLE}\\{A_{i}-E(A_{i}\mid
L_{i})\\}\right]\right)^{2}.$

To estimate the nuisance parameter $E\\{Q_{\tau}(Y\mid A,L)\mid L\\}$, one
would first need to estimate $Q_{\tau}(Y\mid A,L)$. Then, based on these
estimates, one can estimate $E\\{Q_{\tau}(Y\mid A,L)\mid L\\}$ by taking the
quantile estimates as the outcome for a model that predicts its conditional
expectation given $L$. However, when the exposure $A$ is binary, this
expectation reduces to

$E\\{Q_{\tau}(Y\mid A,L)\mid L\\}=Q_{\tau}(Y\mid A=1,L)E(A\mid
L)+Q_{\tau}(Y\mid A=0,L)\\{1-E(A\mid L)\\}.$ (15)

This allows us to estimate this nuisance parameter without additionally
modeling the conditional expectation of $Q_{\tau}(Y\mid A,L)$ given $L$.
Indeed, $Q_{\tau}(Y\mid A,L)$ and $E(A\mid L)$ are nuisance parameters
themselves so they need to be modeled anyway. These models can then be used to
estimate $Q_{\tau}(Y\mid A=1,L)$, $Q_{\tau}(Y\mid A=0,L)$ and $E(A\mid L)$.
Furthermore, targeted learning with cross-fitting can easily be applied in
this case, as

$\hat{E}\\{\widetilde{Q}_{\tau}(Y\mid A,L)\mid L\\}=\widetilde{Q}_{\tau}(Y\mid
A=1,L)\hat{E}(A\mid L)+\widetilde{Q}_{\tau}(Y\mid A=0,L)\\{1-\hat{E}(A\mid
L)\\}.$

On the contrary, decomposition (15) cannot be made for a continuous exposure.
It is then necessary to have an additional model for $E\\{Q_{\tau}(Y\mid
A,L)\mid L\\}$. However, since the retargeting of $\hat{Q}_{\tau}(Y\mid A,L)$
is done on the full data set, it is no longer possible to do proper cross-
fitting for $E\\{Q_{\tau}(Y\mid A,L)\mid L\\}$. In that case, a model to
estimate the conditional expectation of $\widetilde{Q}_{\tau}(Y\mid A,L)$
given $L$ should be trained using only 80% of the data, but
$\widetilde{Q}_{\tau}(Y\mid A,L)$ would already have been targeted using all
data. This makes it no longer possible to train a model for
$E\\{\widetilde{Q}_{\tau}(Y\mid A,L)\mid L\\}$ without involving the remaining
20% of the data. One way to resolve this problem is to do the targeting in
each fold separately. However, forcing equation (11) to hold in every fold
does not guarantee that it holds for the entire data set, especially in
smaller samples. Alternatively, one could choose not to target the quantile
estimator in $E\\{Q_{\tau}(Y\mid A,L)\mid L\\}$. However, instead of repeating
the above targeting algorithm until convergence, we propose to do step 2 only
once, which is justified by the sketched proof in Appendix B. Thus, the
resulting retargeted estimator for $Q_{\tau}(Y\mid A,L)$ will be

$\widetilde{Q}_{\tau}(Y\mid A,L)=\hat{Q}_{\tau}(Y\mid
A,L)+\hat{\epsilon}w(\hat{Q}_{\tau,i}),$

with

$\hat{\epsilon}=\operatorname*{arg\,min}_{\epsilon}\left\lvert\frac{1}{n}\sum_{i=1}^{n}w\big{(}\hat{Q}_{\tau,i}\big{)}\left[\tau-I\left\\{Y_{i}\leq\hat{Q}_{\tau,i}+\epsilon
w\big{(}\hat{Q}_{\tau,i}\big{)}\right\\}\right]\right\rvert.$ (16)

The nuisance parameter $E\\{Q_{\tau}(Y\mid A,L)\mid L\\}$ can then be
estimated based on the following identity:

$\displaystyle\begin{split}E\\{\widetilde{Q}_{\tau}(Y\mid A,L)\mid
L\\}&=E\\{\hat{Q}_{\tau}(Y\mid A,L)+\hat{\epsilon}w(\hat{Q}_{\tau,i})\mid
L\\}\\\ &=E\\{\hat{Q}_{\tau}(Y\mid A,L)\mid
L\\}+\hat{\epsilon}E\left[\frac{A-\hat{E}(A\mid
L)}{\hat{f}\\{\hat{Q}_{\tau}(Y\mid A,L)\mid A,L\\}}\mid L\right].\end{split}$
(17)

Therefore, an estimator can be constructed by estimating
$E\\{\hat{Q}_{\tau}(Y\mid A,L)\mid L\\}$ and $E[\\{A-\hat{E}(A\mid
L)\\}/\hat{f}\\{Q_{\tau}(Y\mid A,L)\mid A,L\\}\mid L]$ separately, possibly
via cross-fitting. In this way, the value of $\epsilon$ can be learned via
(16) after estimating the 2 models in (17) that determine $E\\{Q_{\tau}(Y\mid
A,L)\mid L\\}$, which allows us to do the retargeting on the entire data set.

###  2.4 Parametric quantile regression with variable selection

When using quantile regression, possibly with (stepwise) variable selection,
to estimate the nuisance parameter $Q_{\tau}(Y\mid A,L)$, the debiased
estimator (8) and TMLE estimator (13) can be further simplified. In
particular, assume that $Q_{\tau}(Y\mid A,L)=\beta_{\tau}A+\alpha^{T}L^{*}$,
with $L^{*}$ a subset of variables of $L$. This simplifies the efficient
influence function of $\Psi_{\tau}$ to

$\displaystyle\frac{A-E(A\mid L)}{E[\\{A-E(A\mid
L)\\}^{2}]}\left[\beta\\{A-E(A\mid L)\\}+\frac{\tau-I\\{Y\leq Q_{\tau}(Y\mid
A,L)\\}}{f\\{Q_{\tau}(Y\mid A,L)\mid A,L\\}}-\Psi_{\tau}\\{A-E(A\mid
L)\\}\right].$

The value of $\Psi_{\tau}$ that makes the sample average of the influence
functions zero is then

$\displaystyle\begin{split}\hat{\Psi}_{\tau}&=\hat{\beta}_{\tau}+\frac{1}{n}\sum_{i=1}^{n}\frac{A_{i}-\hat{E}(A_{i}\mid
L_{i})}{\frac{1}{n}\sum_{i=1}^{n}\\{A_{i}-\hat{E}(A_{i}\mid
L_{i})\\}^{2}}\left[\frac{\tau-I\\{Y_{i}\leq\hat{Q}_{\tau}(Y_{i}\mid
A_{i},L_{i})\\}}{\hat{f}\\{\hat{Q}_{\tau}(Y_{i}\mid A_{i},L_{i})\mid
A_{i},L_{i}\\}}\right].\end{split}$ (18)

Here, $\hat{\beta}_{\tau}$ is the estimate of $\beta_{\tau}$ obtained by
fitting a quantile regression model of $Y$ on $A$ and (a subset of) $L$ via
standard software (e.g. the R package ‘quantreg’), possibly with variable
selection. The second term eliminates bias due to variable selection mistakes,
which is needed to achieve valid post-selection inference.

An asymptotically equivalent (one-step) TMLE estimator is

$\displaystyle\begin{split}\hat{\Psi}_{\tau}^{TMLE}&=\hat{\beta}+\frac{1}{n}\sum_{i=1}^{n}\frac{A_{i}-\hat{E}(A_{i}\mid
L_{i})}{\frac{1}{n}\sum_{i=1}^{n}\\{A_{i}-\hat{E}(A_{i}\mid
L_{i})\\}^{2}}\left[\left(\frac{A_{i}-\hat{E}(A_{i}\mid
L_{i})}{\hat{f}\\{\widetilde{Q}_{\tau}(Y_{i}\mid A_{i},L_{i})\mid
A_{i},L_{i}\\}}\right.\right.\\\
&\left.\left.\qquad-\hat{E}\left[\frac{A_{i}-\hat{E}(A_{i}\mid
L_{i})}{\hat{f}\\{\widetilde{Q}_{\tau}(Y_{i}\mid A_{i},L_{i})\mid
A_{i},L_{i}\\}}\mid
L_{i}\right]\right)\hat{\epsilon}+\frac{\tau-I\\{Y_{i}\leq\widetilde{Q}_{\tau}(Y_{i}\mid
A_{i},L_{i})\\}}{\hat{f}\\{\widetilde{Q}_{\tau}(Y_{i}\mid A_{i},L_{i})\mid
A_{i},L_{i}\\}}\right],\end{split}$ (19)

where the last term is approximately zero. More details on the construction of
this estimator can be found in Appendix C.

##  3 What if model (1) is wrong?

Model (1) assumes a linear association between $A$ and the quantiles of $Y$
conditional on $L$, and the absence of $A$-$L$ interactions. Such assumptions
are often made for parsimony and better interpretability, for example, to be
able to summarize the association of interest into one single number. The
downside is that the model may well fail to represent the underlying data-
generating mechanism. Instead of just fitting this model and drawing inference
for $\beta_{\tau}$ as one would do for a model parameter in standard quantile
regression, we have started our analysis by defining a model-free estimand
which reduces to $\beta_{\tau}$ if model (1) is correct, but that still
captures the association of interest if it is misspecified. The latter would,
for instance, be satisfied when the estimand equals some $L$-dependent
weighted average of the estimand which one would report for a subset of
individuals with given $L$ (e.g. in the case of a dichotomous exposure, the
difference in $\tau$-quantile of the outcome between individuals with $A=1$
versus $A=0$ and the same $L$). Building on Vansteelandt and Dukes (2022), we
will focus on the estimand (2) or its equivalent representation (3). For a
dichotomous exposure that is independent of $L$ (e.g. when $A$ is a binary
randomized treatment) this estimand simplifies further to:

$\displaystyle\Psi_{\tau}=E\left\\{Q_{\tau}(Y\mid A=1,L)-Q_{\tau}(Y\mid
A=0,L)\right\\},$ (20)

which equals an expected difference in $\tau$-quantile of the outcome between
individuals with $A=1$ versus $A=0$ and the same level of $L$. The proposed
estimators thus remain interpretable in terms of average differences between
quantiles, even when model (1) is misspecified. This is in stark contrast to
(partially linear) quantile regression that does not provide such a guarantee.

With interest in a causal relationship, one can make use of counterfactual
outcomes $Y^{a}$ to denote the outcome that would have been seen, if the
individual would have had exposure value $A=a$. Under consistency (i.e.
$Y^{a}=Y$) and exchangeability (i.e. when $A$ is independent of the
counterfactual outcome $Y^{a}$ to exposure level $a$,
$Y^{a}\perp\\!\\!\\!\perp A\mid L$), the estimand (20) reduces to

$\displaystyle\Psi_{\tau}=E\left\\{Q_{\tau}(Y^{1}\mid L)-Q_{\tau}(Y^{0}\mid
L)\right\\},$

which is an expected difference in $\tau$-quantile of the counterfactual
outcome $Y^{1}$ versus $Y^{0}$ conditional on $L$.

If $A$ is dichotomous but not randomized, then the estimand reduces to

$\displaystyle\Psi_{\tau}=\frac{E[\pi(L)\\{1-\pi(L)\\}\left\\{Q_{\tau}(Y\mid
A=1,L)-Q_{\tau}(Y\mid A=0,L)\right\\}]}{E[\pi(L)\\{1-\pi(L)\\}]},$ (21)

with $\pi(L)=\text{pr}(A=1\mid L)$, the propensity score. This is a weighted
average of the difference in $\tau$-quantile of the outcome between
individuals with $A=1$ versus $A=0$ and the same level of $L$, with weights
equal to

$\frac{\pi(L)\\{1-\pi(L)\\}}{E[\pi(L)\\{1-\pi(L)\\}]}.$

This choice of weights is motivated by settings where interest lies in the
causal effect of $A$. One then wishes to upweight those strata where ${\rm
Var}(A\mid L)=\pi(L)\\{1-\pi(L)\\}$ is relatively large. This corresponds to
giving the highest weight to covariate regions where both exposed and
unexposed subjects are found. Moreover, in strata of $L$ where all subjects
are exposed (or unexposed), those subjects will not contribute to
$\Psi_{\tau}$. Those strata do not contain pertinent information about the
causal effect of $A$ and will be given a weight of zero.

As above, under consistency and conditional exchangeability, (21) reduces to a
weighted average of the difference in $\tau$-quantile of $Y^{1}$ versus
$Y^{0}$ conditional on $L$:

$\displaystyle\Psi_{\tau}=\frac{E\left[\pi(L)\\{1-\pi(L)\\}\left\\{Q_{\tau}(Y^{1}\mid
L)-Q_{\tau}(Y^{0}\mid L)\right\\}\right]}{E[\pi(L)\\{1-\pi(L)\\}]}.$

For exposures that are not necessarily binary, this estimand can be
generalized by (2), which is constructed as follows. Consider drawing a sample
of individuals, and for each, randomly drawing an individual with the same
$L$. A population least squares regression of the difference between their
conditional quantiles on their difference in exposure values results in (2),
or equivalently, in (3). This is closely related to the idea of matching
(Rosenbaum and Rubin, 1983; Rosenbaum, 1989) based on $L$.

Further insight can be obtained by considering the case where

$Q_{\tau}(Y\mid A,L)=\beta_{\tau}(L)A+\omega_{\tau}(L),$

where the coefficient $\beta_{\tau}$ is now allowed to be $L$-dependent. By
substituting the right-hand side in (3), the estimand can easily be shown to
reduce to

$\Psi_{\tau}=\frac{E\left\\{{\rm Var}(A\mid
L)\beta_{\tau}(L)\right\\}}{E\left\\{{\rm Var}(A\mid L)\right\\}}.$

This is a weighted average of the $L$-dependent coefficients
$\beta_{\tau}(L)$, which thus captures the conditional association of
interest.

##  4 Simulation experiments

To investigate the performance of the proposed estimators, we use four
simulation experiments. In all experiments, we report Monte Carlo bias and
standard deviation (SD), as well as the averaged estimated standard error (SE)
and coverage of corresponding 95% Wald confidence intervals. For all
experiments, we show results for 1000 simulation runs for the oracle estimator
(Oracle), obtained by fitting a correctly specified standard quantile
regression model that includes the exposure $A$ and all necessary confounders,
the naive plug-in (Plug-in) estimator (4), the DML estimator (8) based on
quantile random forests without cross-fitting (DML) and with 5-fold cross-
fitting (DML-CF), and the targeted learning estimator (13) without cross-
fitting (TMLE) and with 5-fold cross-fitting (TMLE-CF). The sample size is
fixed at $n=500$ unless otherwise stated. The standard error of the plug-in
estimator is always calculated similarly as for the debiased machine learning
estimator, but removing the last term of (8) from the efficient influence
function. More experiment-specific details are provided below.

###  4.1 Experiment 1: Homoscedasticity, heteroscedasticity, and continuity

This experiment consisted of 3 simulation settings, where the data generating
mechanism was each time slightly changed. We studied the proposed estimators’
behaviour for a binary exposure in the first simulation setting, inspired by a
setting used in Sun et al. (2021). Some minor adjustments were made to prevent
the coefficient of determination ($R^{2}$) from exceeding 0.8, as this is
often unrealistic. In setting 1, we considered a homoscedastic outcome where
the exposure effect had the same magnitude for all quantiles. Thus, there was
a distribution shift between exposed and unexposed individuals. In setting 2,
we adjusted the setting slightly so that the outcome variability differed for
both exposure groups. This implies that the exposure effect is different for
different values of $\tau$. In setting 3, the exposure was continuous. The
results are shown in Table 1.

In the first setting, the data was generated as follows. Four confounding
variables were generated according to

$\displaystyle\begin{pmatrix}L_{1}\\\ L_{2}\\\ L_{3}\\\
L_{4}\end{pmatrix}\sim\mathcal{N}\left\\{\begin{pmatrix}0\\\ 0\\\ 0\\\
0\end{pmatrix},\begin{pmatrix}1&0.5&0.2&0.3\\\ 0.5&1&0.7&0\\\ 0.2&0.7&1&0\\\
0.3&0&0&1\end{pmatrix}\right\\}.$ (22)

The dichotomous exposure $A$ follows a Bernoulli distribution conditional on
$L=(L_{1},L_{2},L_{3},L_{4})$, with true propensity score given by
$\text{pr}(A=1\mid L)=\text{expit}(-0.5+0.2L_{1}-0.4L_{2}-0.4L_{3}+0.2L_{4})$.

In the homoscedastic setting, the outcome $Y$ was generated as

$\displaystyle Y=1+A+\sin(L_{1})+L_{2}^{2}+L_{3}+L_{4}+L_{3}L_{4}+\epsilon,$
(23)

where $\epsilon$ follows a Gamma($k,\theta$) distribution with shape parameter
$k=1$ and scale parameter $\theta=2$ and is independent from $A$ and $L$.
Under heteroscedasticity, we chose $\epsilon\sim\text{Gamma}(1,2+A)$.

For all estimators except Oracle, quantile random forests (via the ‘grf’
package described in Athey et al. (2019)) were used to learn $Q_{\tau}(Y\mid
A,L)$. To learn $E(A\mid L)$, SuperLearner was used; this is an ensemble
learner described in van der Laan et al. (2007). The SuperLearner library
included generalized linear models with and without elastic net
regularization, generalized additive models, and random forests. These
estimators were also used to estimate $E\\{Q_{\tau}(Y\mid A,L)\mid L\\}$ based
on (15). To obtain the nuisance parameter $f\\{Q_{\tau}(Y\mid A,L)\mid
A,L\\}$, we estimated the conditional density of $Y$ given $A$ and $L$ via the
‘fk_density’ function from the R package ‘FKSUM’ (Hofmeyr, 2021, 2022). This
function gives a fast approximation of a univariate density based on recursive
computations in evaluation points of choice. Since our interest lies in
estimating a conditional density function, we estimated the density of
$Y-Q_{\tau}(Y\mid A,L)$ and evaluated it at 0. This estimation procedure
assumes that the distribution of $Y-Q_{\tau}(Y\mid A,L)$ is independent of $A$
and $L$.

|  | $\tau=0.5$ |  | $\tau=0.75$ |  | $\tau=0.9$
---|---|---|---|---|---|---
Setting | estimator | bias | SD | SE | Cov |  | bias | SD | SE | Cov |  | bias | SD | SE | Cov
| Oracle | -0.17 | 19 | 20 | 96.6 |  | -0.51 | 33 | 35 | 96.3 |  | -1.1 | 56 | 60 | 96.0
| Plug-in | -70 | 12 | 1.5 | 0.1 |  | -72 | 14 | 1.8 | 0.5 |  | -64 | 22 | 3.6 | 1.6
| DML | -15 | 19 | 19 | 85.9 |  | -35 | 29 | 25 | 63.7 |  | -56 | 32 | 22 | 33.0
| DML-CF | -16 | 18 | 18 | 84.2 |  | -24 | 28 | 25 | 79.2 |  | -31 | 43 | 35 | 77.4
| TMLE | -41 | 21 | 19 | 45.2 |  | -71 | 32 | 26 | 27.0 |  | -110 | 28 | 45 | 12.6
1 | TMLE-CF | 1.2 | 22 | 25 | 97.2 |  | 2.8 | 39 | 37 | 93.5 |  | 14 | 68 | 63 | 91.4
| Oracle | -0.59 | 26 | 26 | 93.7 |  | -2.7 | 45 | 46 | 93.6 |  | -6.8 | 75 | 78 | 94.4
| Plug-in | -89 | 26 | 2.7 | 0.6 |  | -120 | 42 | 4.1 | 1.1 |  | -160 | 60 | 6.8 | 1.5
| DML | -15 | 27 | 22 | 83.0 |  | -45 | 46 | 30 | 57.8 |  | -120 | 66 | 28 | 16.5
| DML-CF | -20 | 26 | 20 | 75.6 |  | -43 | 43 | 30 | 61.7 |  | -78 | 69 | 42 | 49.6
| TMLE | -45 | 27 | 23 | 49.4 |  | -76 | 49 | 32 | 38.0 |  | -195 | 67 | 34 | 4.0
2 | TMLE-CF | 5.0 | 29 | 28 | 93.9 |  | 10 | 50 | 42 | 88.7 |  | 13 | 88 | 70 | 87.6
| Oracle | -0.13 | 3.5 | 3.6 | 95.6 |  | -0.29 | 5.8 | 6.3 | 95.8 |  | 0.10 | 10 | 11 | 94.6
| Plug-in | -17 | 6.4 | 1.6 | 0.5 |  | -25 | 7.9 | 1.9 | 0.2 |  | -39 | 11 | 2.1 | 0.0
| DML | -7.9 | 6.8 | 3.6 | 45.1 |  | -16 | 8.5 | 4.6 | 21.6 |  | -34 | 11 | 4.3 | 1.7
| DML-CF | -6.2 | 3.9 | 3.1 | 48.5 |  | -10 | 6.0 | 4.4 | 40.0 |  | -19 | 9.1 | 6.2 | 23.0
| TMLE | -11 | 7.1 | 3.6 | 29.3 |  | -18 | 9.3 | 4.9 | 18.9 |  | -41 | 12 | 5.2 | 1.6
3 | TMLE-CF | -1.1 | 4.4 | 4.2 | 92.9 |  | 0.016 | 7.7 | 6.3 | 90.0 |  | 1.2 | 14 | 10 | 85.3
Table 1: Simulation results for experiment 1: sample size $n=500$, quantile
$\tau$, Monte Carlo bias, Monte Carlo standard deviation (SD), average of the
influence function based standard errors (SE) and coverage of 95% Wald
confidence intervals (Cov). All values have been multiplied by $10^{2}$.

Estimating the conditional density could alternatively be done via the
‘haldensify’ package (Hejazi et al., 2022a, b, c), which uses the highly
adaptive lasso (HAL) (Benkeser and Van Der Laan, 2016; van der Laan, 2017).
This conditional density could then be evaluated in the estimated quantile
without making any additional assumptions. However, this method is much more
time-consuming and did not outperform the simpler estimation method. The
corresponding results can be found in Appendix D.

In the third setting of this experiment, we studied the proposed estimators’
finite sample behaviour for a continuous exposure. Here, $A$ was generated
according to a normal distribution $\mathcal{N}(\mu,2^{2})$, with
$\mu=-0.5+L_{1}-2L_{2}-2L_{3}+L_{4}$. The outcome $Y$ was generated from (23),
where $\epsilon\sim\text{Gamma}(1,4)$ is independent from $A$ and $L$. The
nuisance parameters were all estimated as before, except for
$E\\{Q_{\tau}(Y\mid A,L)\mid L\\}$ as the decomposition (15) is no longer
possible for continuous $A$. To estimate this parameter, we again used
SuperLearner with the same library as before, using the random forest estimate
of $Q_{\tau}(Y\mid A,L)$ as outcome. For the TMLE-CF estimator, the nuisance
parameters described at the end of Section 2.3 are also estimated using
SuperLearner.

The simulations all demonstrate the poor performance of the naive plug-in
estimator. It has a high bias and the standard error is poorly estimated and
overly optimistic, leading to a coverage of almost 0. When comparing the
debiased machine learning and targeted learning estimators with versus without
cross-fitting, the advantage of cross-fitting becomes clearly visible. In
general, the estimators with cross-fitting perform better, especially for
$\tau=0.9$ and in the setting with continuous exposure. The use of cross-
fitting leads to a bias reduction and to estimated standard errors that better
approximate the Monte Carlo standard deviation, resulting in a better
coverage. For $\tau=0.5$, we sometimes see the estimators without cross-
fitting performing slightly better in terms of bias, likely because the
nuisance parameters are estimated on larger sample sizes. When comparing DML-
CF with TMLE-CF, we often observe the latter having a larger SD. However, in
terms of bias and coverage, we always observe targeted learning outperforming
the other estimator. As expected, all estimators perform worse for more
extreme quantiles.

###  4.2 Experiment 2: Extreme propensity scores

In this experiment, we generated a setting where the propensity scores could
take extreme values (i.e. values close to 0 or 1), leading to little overlap
between exposed and unexposed individuals. Therefore, we generated a
dichotomous exposure $A$, with $\text{pr}(A=1\mid
L)=\text{expit}(-0.5+0.2L_{1}-0.4L_{2}-0.4L_{3}+0.2L_{4}+0.5L_{1}^{2}-0.5L_{2}^{2}+0.5L_{3}L_{4})$,
where $L$ follows the same multivariate normal distribution (22) as in
experiment 1. To gain insight into these propensity scores, histograms are
provided in Appendix D. The outcome was generated according to (23), with
$\epsilon\sim\text{Gamma}(1,3)$. We evaluated the same estimators as in
experiment 1 and included a naive (misspecified) quantile regression model
(QR) that only included main effects. The results are shown in Table 2.

Sample |  | $\tau=0.5$ |  | $\tau=0.75$ |  | $\tau=0.9$
---|---|---|---|---|---|---
size | estimator | bias | SD | SE | Cov |  | bias | SD | SE | Cov |  | bias | SD | SE | Cov
| Oracle | -0.30 | 29 | 30 | 94.6 |  | 0.68 | 51 | 52 | 94.8 |  | -3.0 | 92 | 91 | 94.2
| QR | 16 | 36 | 36 | 92.3 |  | 33 | 55 | 58 | 92.7 |  | 44 | 91 | 95 | 94.1
| Plug-in | -72 | 14 | 1.5 | 0.0 |  | -74 | 17 | 2.0 | 0.4 |  | -62 | 36 | 4.7 | 5.8
| DML | -12 | 32 | 33 | 93.1 |  | -37 | 44 | 41 | 82.0 |  | -56 | 54 | 35 | 54.0
| DML-CF | -21 | 27 | 25 | 85.3 |  | -35 | 40 | 36 | 80.5 |  | -39 | 63 | 46 | 74.7
| TMLE | -59 | 33 | 33 | 55.0 |  | -97 | 45 | 43 | 35.1 |  | -110 | 97 | 45 | 31.8
$500$ | TMLE-CF | -6.5 | 32 | 34 | 95.4 |  | -13 | 55 | 52 | 92.9 |  | -14 | 100 | 87 | 91.2
| Oracle | -0.089 | 20 | 21 | 95.6 |  | 0.37 | 36 | 37 | 95.0 |  | -3.0 | 61 | 63 | 95.3
| QR | 17 | 25 | 26 | 91.1 |  | 33 | 40 | 40 | 86.4 |  | 43 | 64 | 65 | 90.8
| Plug-in | -66 | 12 | 1.1 | 0.0 |  | -69 | 16 | 1.5 | 0.2 |  | -60 | 27 | 3.2 | 2.1
| DML | -14 | 21 | 22 | 90.5 |  | -34 | 32 | 29 | 74.3 |  | -57 | 36 | 23 | 34.0
| DML-CF | -22 | 19 | 17 | 73.6 |  | -32 | 29 | 24 | 68.1 |  | -41 | 43 | 29 | 59.7
| TMLE | -59 | 22 | 22 | 24.1 |  | -100 | 33 | 30 | 9.7 |  | -130 | 86 | 29 | 17.2
$1000$ | TMLE-CF | -6.8 | 22 | 23 | 94.3 |  | -9.5 | 39 | 37 | 93.0 |  | -14 | 69 | 62 | 90.7
Table 2: Simulation results for experiment 2: sample size $n$, quantile
$\tau$, Monte Carlo bias, Monte Carlo standard deviation (SD), average of the
influence function based standard errors (SE) and coverage of 95% Wald
confidence intervals (Cov). All values have been multiplied by $10^{2}$.

The proposed estimators show similar behaviour in this experiment as in the
first experiment. The best performance is again achieved by the TMLE-CF
estimator. It clearly performs better than a naive quantile regression,
particularly in terms of bias. When increasing the sample size from 500 to
1000, TMLE-CF does not significantly change performance. However, QR’s bias
does not decrease, unlike its standard error, leading to a decrease in
coverage of 95% confidence intervals. The TMLE estimator also has a large
bias, which is likely caused by overfitting since the bias of TMLE-CF is much
smaller. As in the first experiment, cross-fitting provides a major
improvement.

###  4.3 Experiment 3: Randomized dichotomous exposure

In this experiment, we first chose the exposure to be randomized with
$\text{pr}(A=1)=0.5$, independently from all other variables. Confounders
follow the same multivariate normal distribution (22) as in experiment 1. The
outcome was generated from (23), again with $\epsilon\sim\text{Gamma}(1,2)$.
We compared our estimators with 2 estimators recently proposed in Athey et al.
(2023). Both estimators focus on comparing marginal quantiles
$Q_{\tau}(Y^{a})$ of counterfactual distributions $Y^{a}$, $a=0,1$. They are
designed for a constant quantile treatment effect between 2 groups, meaning
that there is a location shift between the outcome distribution of exposed
versus unexposed subjects. Athey et al. (2023) estimate so-called weighted
average quantile treatment effects

$\displaystyle\int_{0}^{1}\left\\{Q_{\tau}(Y^{1})-Q_{\tau}(Y^{0})\right\\}\,\mathrm{d}W(\tau),$
(24)

for certain weights $W(\tau)$ that integrate to one, and $W(0)=0$, $W(1)=1$.
The first estimator (IFB) is an influence function based estimator, derived
from the efficient influence function of (24). It is a semiparametric
efficient estimator for (24). The second estimator (WAQ) is a weighted average
quantile estimator. It is a plug-in estimator for (24):

$\displaystyle\int_{0}^{1}\left\\{\hat{Q}_{\tau}(Y^{1})-\hat{Q}_{\tau}(Y^{0})\right\\}w_{f}(\tau)\,\mathrm{d}\tau,$

where the weights $w_{f}(\tau)$ are such that they provide an efficient
L-estimator under the location shift assumption and the assumption that the
density $f$ of the counterfactual outcome is known (see Athey et al. (2023)
for details). Since these estimators assume a location shift, they are only
implemented for estimating the median. The results are shown in Table 3
(left).

For the proposed estimators we see similar results as in experiment 1.
Targeted learning with cross-fitting clearly outperforms the other estimators
with a bias of 0.013 and coverage of 95.1%. Despite our proposal being
assumption-lean, it still performs slightly better than the estimators defined
in Athey et al. (2023), despite all assumptions being met in this simulation
experiment.

estimator | bias | SD | SE | Cov |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---|---
Oracle | 0.0071 | 0.18 | 0.19 | 95.5 |  |  |  |  |  |  |
Plug-in | -0.70 | 0.12 | 0.011 | 0.0 |  |  |  |  |  |  |
DML | -0.16 | 0.20 | 0.18 | 80.3 |  | estimator | bias | SD | SE | Cov |
DML-CF | -0.18 | 0.19 | 0.17 | 76.1 |  | TMLE-CF | 0.050 | 0.29 | 0.28 | 93.9 |
TMLE | -0.42 | 0.21 | 0.18 | 37.6 |  | IFB | -0.21 | 0.23 | 0.25 | 88.0 |
TMLE-CF | 0.013 | 0.23 | 0.23 | 95.1 |  | WAQ | -0.19 | 0.25 | 0.25 | 88.5 |
IFB | 0.016 | 0.22 | 0.23 | 96.4 |  |  |  |  |  |  |
WAQ | 0.018 | 0.22 | 0.23 | 95.8 |  |  |  |  |  |  |
Table 3: Simulation results for (left) randomized dichotomous exposure and
(right) setting 2 of experiment 1: sample size $n=500$, quantile $\tau=0.5$,
Monte Carlo bias, Monte Carlo standard deviation (SD), average of the
influence function based standard errors (SE) and coverage of 95% Wald
confidence intervals (Cov).

However, when applying the two estimators from Athey et al. (2023) in the
heteroscedasticity setting from experiment 1, we observe a much worse
performance (see Table 3 (right)). Here the location shift assumption is
violated, leading to an inflated bias. As expected, our proposal outperforms
their estimators. The results of the TMLE-CF estimator are taken from Table 1
and included in Table 3 (right) for comparison.

###  4.4 Experiment 4: Variable selection in higher dimensions

In this simulation experiment, we considered a covariate vector of dimension
50, based on the data-generating mechanism in Belloni et al. (2013a). In
particular, we generated a 50-dimensional covariate vector
$L\sim\mathcal{N}(0,\Sigma)$, where $\Sigma$ was an autoregressive correlation
matrix with correlation parameter 0.5. A continuous exposure was generated
according to a normal distribution with mean $\sum_{k=1}^{10}L_{k}/k$ and unit
residual variance. The outcome $Y$ was generated from
$Y=A+\sum_{k=1}^{5}L_{k}/k+\sum_{k=11}^{15}L_{k}/(k-10)+\epsilon$, with
$\epsilon\sim\mathcal{N}(0,2^{2})$ independent from $A$ and $L$. Here, we also
compare sample size $n=250$ with $n=500$.

In this experiment, we compare the performance of a quantile regression with
AIC-based backward stepwise variable selection (QRvs) with the TMLE estimator
defined in (19) with 5-fold cross-fitting (TMLEvs-CF). The results are shown
in Table 4.

Appendix D contains the entire table with additional results for the
estimators used in experiment 1, the DML estimator defined in (18) where
quantile random forests are replaced by quantile regression with AIC-based
backward stepwise variable selection and the TMLE estimator defined in (19)
without cross-fitting.

|  |  | $\mathbf{n=250}$ |  | $\mathbf{n=500}$
---|---|---|---|---|---
Quantile | estimator |  | bias | SD | SE | Cov |  | bias | SD | SE | Cov
$\tau=0.5$ | Oracle |  | 0.66 | 15 | 16 | 96.0 |  | 0.25 | 10 | 11 | 95.7
QRvs |  | 2.4 | 17 | 16 | 91.9 |  | 0.36 | 12 | 11 | 92.8
TMLEvs-CF |  | 1.3 | 20 | 19 | 94.3 |  | -0.32 | 13 | 12 | 94.8
$\tau=0.75$ | Oracle |  | -0.016 | 16 | 18 | 96.7 |  | -0.54 | 12 | 12 | 95.0
QRvs |  | -1.4 | 19 | 17 | 91.3 |  | 0.0095 | 13 | 12 | 93.2
TMLEvs-CF |  | 0.52 | 21 | 21 | 94.0 |  | -0.28 | 14 | 13 | 93.5
$\tau=0.9$ | Oracle |  | -0.88 | 21 | 22 | 95.7 |  | 0.39 | 14 | 15 | 95.8
QRvs |  | 0.065 | 23 | 20 | 89.7 |  | 0.094 | 16 | 14 | 90.8
TMLEvs-CF |  | 0.21 | 24 | 24 | 94.4 |  | -0.18 | 16 | 16 | 93.7
Table 4: Simulation results for variable selection: sample size $n$, quantile
$\tau$, Monte Carlo bias, Monte Carlo standard deviation (SD), average of the
influence function based standard errors (SE) and coverage of 95% Wald
confidence intervals (Cov). All values have been multiplied by $10^{2}$.

We observe that the standard quantile regression with variable selection
(QRvs) performs well in terms of bias. However, there is a slight
underestimation of the standard errors, leading to a coverage slightly below
95%. On the contrary, the TMLE estimator achieves a coverage of approximately
94%. As expected, its standard errors are more accurate because it takes the
excess uncertainty of the variable selection procedure into account.

##  5 Data analysis

To illustrate our proposal, we reanalyze an observational study to estimate
annual health care and lost productivity costs associated with excess weight
among the adult population in Belgium, using national health data. Gorasso et
al. (2022) estimated the average difference in health care costs between
individuals with underweight (BMI $<18.5$) and normal weight ($18.5\leq$ BMI
$<25$) versus overweight ($25\leq$ BMI $<30$) and obesity (BMI $\geq 30$)
using a linear main effect model. They used a double-selection approach to
select potential confounders (Belloni et al., 2013b) and estimated a
population level effect using g-computation (Robins, 1986). Candidate
variables included age groups, gender, household educational level, household
level of income, and some behavioral risk factors such as smoking, alcohol
misuse, unhealthy eating behavior, and physical inactivity.

Since health care costs are right-skewed, the average effect of overweight and
obesity could be highly influenced by extreme costs. Therefore, analyzing
quantiles could add more insight. We will use the proposed TMLE estimator (19)
with 5-fold cross-fitting to assess the conditional association between health
care cost and body weight (dichotomized as underweight and normal weight vs.
overweight and obesity), considering the same set of candidate baseline
variables as Gorasso et al. (2022). This estimator makes use of parametric
quantile regression with AIC-based backward stepwise variable selection. We
will compare the results with a traditional parametric quantile model that
uses the same variable selection procedure. All analyses include sampling
weights to adjust for the multi-stage sampling design (Demarest et al., 2013).

Figure 1: Effect estimate of overweight and obesity on annual health care
costs per capita, with 95% confidence intervals for the proposed method (ALQR)
and a traditional quantile regression approach (Traditional QR).

With access to complete data of 3925 individuals, the traditional parametric
quantile regression with stepwise variable selection delivers an estimate of
the median incremental health care cost per capita of a high BMI of €187.70
(SE: 53.95) and a 90th percentile equal to €521.13 (SE: 259.05). The proposed
TMLE approach delivers an estimated median incremental health care cost per
capita of excess body weight of €205.09 (SE: 75.94) and an estimate of
€1142.42 (SE: 220.15) for the 90th percentile. The fact that our proposal
yields different results may be caused by misspecification of the parametric
quantile regression model, which did not include interactions involving BMI
class. For instance, there is a significant interaction between BMI class and
education level ($p=0.0023$). Results for more quantiles are shown in Fig. 1.
We observe that the difference in annual health care cost for high versus low
body weight is much larger when looking at higher quantiles of the conditional
outcome distribution. A sensitivity analysis, which examines additional
variability for extreme quantiles related to randomness in the cross-fitting,
is included in Appendix E.

##  6 Discussion

In this article, inspired by Vansteelandt and Dukes (2022), we have focused on
providing a scalar quantile-based summary of the conditional association
between an outcome and an exposure, given covariates. Our results are
primarily of interest in settings where averages are less informative and
interest rather lies in median (or other quantile) effects (e.g. for heavy-
tailed or skewed data). The main motivation to focus on conditional quantiles
is because marginal quantiles of counterfactual outcomes often demand strong
extrapolations by considering a hypothetical scenario where everyone got (or
did not get) treatment. This is problematic when some covariate groups are
highly unlikely to get treatment. We accommodate this by letting those
covariate groups contribute little to our estimand via a very small weight, as
is also common in standard regression methods. These groups contain very
little information on the exposure effect and there may often be less interest
in the effect for those individuals. However, we faced challenges because of
the conditional density $f\\{Q_{\tau}(Y\mid A,L)\mid A,L\\}$ appearing as a
nuisance parameter in the estimator. Estimation of this density is generally a
difficult task. Like Zhong and Wang (2023), we assume homoscedastic residuals
to estimate it. The ‘FKSUM’ estimator for instance assumes this conditional
density to be equal for all subjects. One way to sidestep this assumption is
to estimate the (marginal) distribution of $\\{Y-E(Y\mid A,L)\\}{\rm
Var}(Y\mid A,L)^{-1/2}$ and then transform an estimate of this distribution to
estimate $f\\{Q_{\tau}(Y\mid A,L)\mid A,L\\}$. However, this approach resulted
in an extremely unstable behaviour of our estimator for $\Psi_{\tau}$ (not
shown) because of estimates of this density close to zero. The alternative HAL
estimator (Benkeser and Van Der Laan, 2016; van der Laan, 2017) is more
flexible, and does not suffer from instability, but is computationally very
expensive. Since the performance of the TMLE estimator did not significantly
improve when using the latter density estimator (even in the absence of a
location shift), we opted for the simpler one. A broader exploration of
different (conditional) density estimators may lead to further improvements in
the proposed estimator. Unfortunately, the existing literature on inverse
density weighting is not directly applicable (Huling et al., 2023), because we
need to evaluate the density in a quantile instead of in the observed
outcomes. In future work, we will examine if Auto-DML (Chernozhukov et al.,
2020, 2022) could offer a solution to directly estimate the inverse density
rather than first estimating the density and then taking the reciprocal.

From the simulation results for more extreme quantiles, it was often observed
that the average estimated standard error is less than the Monte Carlo
standard deviation of the estimator. This could be the result of ignoring the
randomness in cross-fitting, for which rank-transformed subsampling (Guo and
Shah, 2023) could offer improvements.

Our proposed estimators allow for estimating the nuisance parameters through
data-adaptive techniques. To estimate $Q_{\tau}(Y\mid A,L)$, we utilized
either quantile regression forests or parametric quantile regression with
variable selection. A potential improvement could be achieved by using an
ensemble learner, e.g. by constructing an ensemble through the R package
‘SuperLearner’, including the aforementioned methods in its library,
potentially supplemented with other estimation methods (e.g. Lasso penalized
quantile regression). In this way, the user no longer needs to manually choose
which algorithm to use. However, to our knowledge such ensemble learners for
quantile estimation are currently unavailable.

Finally, it could be useful to generalize the estimand (3) so that it includes
a link function and replace model (1) by $g\\{Q_{\tau}(Y\mid
A,L)\\}=\beta_{\tau}A+\omega_{\tau}(L)$, with link function $g$. More details
about this proposal can be found in Appendix F.

## References

  * Koenker and Bassett [1978] Roger Koenker and Gilbert Bassett.  Regression Quantiles.  _Econometrica_ , 46(1):33–50, 1978.  ISSN 00129682, 14680262.  URL http://www.jstor.org/stable/1913643.
  * Koenker [2005] Roger Koenker.  _Quantile Regression_ , volume 38.  Cambridge University Press, 2005.
  * Koenker [2017] Roger Koenker.  Quantile regression: 40 years on.  _Annual Review of Economics_ , 9(1):155–176, 2017.  doi: 10.1146/annurev-economics-063016-103651.  URL https://doi.org/10.1146/annurev-economics-063016-103651.
  * Koenker et al. [2017] Roger Koenker, Victor Chernozhukov, Xuming He, and Limin Peng.  _Handbook of quantile regression_.  Chapman & Hall/CRC, New York, 1st edition, 2017.
  * Koenker and Hallock [2001] Roger Koenker and Kevin F. Hallock.  Quantile Regression.  _Journal of Economic Perspectives_ , 15(4):143–156, 2001.  doi: 10.1257/jep.15.4.143.  URL https://www.aeaweb.org/articles?id=10.1257/jep.15.4.143.
  * Yu et al. [2003] Keming Yu, Zudi Lu, and Julian Stander.  Quantile regression: applications and current research areas.  _Journal of the Royal Statistical Society: Series D (The Statistician)_ , 52(3):331–350, 2003.  doi: https://doi.org/10.1111/1467-9884.00363.  URL https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/1467-9884.00363.
  * Das et al. [2019] Kiranmoy Das, Martin Krzywinski, and Naomi Altman.  Quantile regression.  _Nature Methods_ , 16(6):451–452, 2019.
  * Hernán and Robins [2020] Miguel A Hernán and James M Robins.  _Causal inference: What If_.  Boca Raton: Chapman & Hall/CRC, 2020.
  * Doksum [1974] Kjell Doksum.  Empirical Probability Plots and Statistical Inference for Nonlinear Models in the Two-Sample Case.  _The Annals of Statistics_ , 2(2):267–277, 1974.  ISSN 00905364.  URL http://www.jstor.org/stable/2958036.
  * Lehmann and D’Abrera [1975] E. L. Lehmann and H. J. D’Abrera.  _Nonparametrics: Statistical methods based on ranks_.  Holden-Day, 1975.
  * Vansteelandt and Dukes [2022] Stijn Vansteelandt and Oliver Dukes.  Assumption-lean Inference for Generalised Linear Model Parameters (with discussion).  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 84(3):657–685, 07 2022.  ISSN 1369-7412.  doi: 10.1111/rssb.12504.  URL https://doi.org/10.1111/rssb.12504.
  * Breiman [2001] Leo Breiman.  Statistical Modeling: The Two Cultures (with comments and a rejoinder by the author).  _Statistical Science_ , 16(3):199–231, 2001.  doi: 10.1214/ss/1009213726.  URL https://doi.org/10.1214/ss/1009213726.
  * Leeb and Pötscher [2005] Hannes Leeb and Benedikt M. Pötscher.  Model selection and inference: Facts and fiction.  _Econometric Theory_ , 21(1):21–59, 2005.  doi: 10.1017/S0266466605050036.
  * Melly [2006] Blaise Melly.  Estimation of counterfactual distributions using quantile regression.  _University of St.Gallen_ , 68, 01 2006.
  * Firpo [2007] Sergio Firpo.  Efficient Semiparametric Estimation of Quantile Treatment Effects.  _Econometrica_ , 75(1):259–276, 2007.  doi: https://doi.org/10.1111/j.1468-0262.2007.00738.x.  URL https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1468-0262.2007.00738.x.
  * Chernozhukov et al. [2013] Victor Chernozhukov, Iván Fernández-Val, and Blaise Melly.  Inference on Counterfactual Distributions.  _Econometrica_ , 81(6):2205–2268, 2013.  doi: https://doi.org/10.3982/ECTA10582.  URL https://onlinelibrary.wiley.com/doi/abs/10.3982/ECTA10582.
  * Donald and Hsu [2014] Stephen G. Donald and Yu-Chin Hsu.  Estimation and inference for distribution functions and quantile functions in treatment effect models.  _Journal of Econometrics_ , 178:383–397, 2014.  ISSN 0304-4076.  doi: https://doi.org/10.1016/j.jeconom.2013.03.010.  URL https://www.sciencedirect.com/science/article/pii/S0304407613001826.
  * Caracciolo and Furno [2017] Francesco Caracciolo and Marilena Furno.  Quantile treatment effect and double robust estimators: an appraisal on the italian labor market.  _Journal of Economic Studies_ , 4(2):585–604, 2017.
  * Athey et al. [2023] Susan Athey, Peter J Bickel, Aiyou Chen, Guido W Imbens, and Michael Pollmann.  Semi-parametric estimation of treatment effects in randomised experiments.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 07 2023.  ISSN 1369-7412.  doi: 10.1093/jrsssb/qkad072.  URL https://doi.org/10.1093/jrsssb/qkad072.
  * Ai et al. [2022] Chunrong Ai, Oliver Linton, and Zheng Zhang.  Estimation and inference for the counterfactual distribution and quantile functions in continuous treatment models.  _Journal of Econometrics_ , 228(1):39–61, 2022.  ISSN 0304-4076.  doi: https://doi.org/10.1016/j.jeconom.2020.12.009.  URL https://www.sciencedirect.com/science/article/pii/S0304407621000543.  Annals Issue: In Honor of Ron Gallant.
  * Díaz [2017] Iván Díaz.  Efficient estimation of quantiles in missing data models.  _Journal of Statistical Planning and Inference_ , 190:39–51, 2017.  ISSN 0378-3758.  doi: https://doi.org/10.1016/j.jspi.2017.05.001.  URL https://www.sciencedirect.com/science/article/pii/S0378375817300782.
  * Belloni et al. [2017] Alexandre Belloni, Victor Chernozhukov, Ivan Fernandez-Val, and Christian Hansen.  Program evaluation and causal inference with high-dimensional data.  _Econometrica_ , 85(1):233–298, 2017.
  * Kallus et al. [2019] Nathan Kallus, Xiaojie Mao, and Masatoshi Uehara.  Localized debiased machine learning: Efficient inference on quantile treatment effects and beyond.  _arXiv preprint arXiv:1912.12945_ , 2019.
  * Chernozhukov et al. [2018] Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins.  Double/debiased machine learning for treatment and structural parameters.  _The Econometrics Journal_ , 21(1):C1–C68, 01 2018.  ISSN 1368-4221.  doi: 10.1111/ectj.12097.  URL https://doi.org/10.1111/ectj.12097.
  * Lee [2003] Sokbae Lee.  Efficient semiparametric estimation of a partially linear quantile regression model.  _Econometric Theory_ , 19(1):1–31, 2003.  doi: 10.1017/S0266466603191013.
  * Sun [2005] Yiguo Sun.  Semiparametric efficient estimation of partially linear quantile regression models.  _Annals of Economics and Finance_ , 6(1):105, 2005.
  * Wu et al. [2010] Tracy Z. Wu, Keming Yu, and Yan Yu.  Single-index quantile regression.  _Journal of Multivariate Analysis_ , 101(7):1607–1621, 2010.  ISSN 0047-259X.  doi: https://doi.org/10.1016/j.jmva.2010.02.003.  URL https://www.sciencedirect.com/science/article/pii/S0047259X10000333.
  * Wu and Yu [2014] Chaojiang Wu and Yan Yu.  Partially linear modeling of conditional quantiles using penalized splines.  _Computational Statistics & Data Analysis_, 77:170–187, 2014.  ISSN 0167-9473.  doi: https://doi.org/10.1016/j.csda.2014.02.020.  URL https://www.sciencedirect.com/science/article/pii/S0167947314000607.
  * Lv et al. [2015] Yazhao Lv, Riquan Zhang, Weihua Zhao, and Jicai Liu.  Quantile regression and variable selection of partial linear single-index model.  _Annals of the Institute of Statistical Mathematics_ , 67(2):375–409, 2015.
  * Sherwood and Wang [2016] Ben Sherwood and Lan Wang.  Partially linear additive quantile regression in ultra-high dimension.  _The Annals of Statistics_ , 44(1):288–317, 2016.  doi: 10.1214/15-AOS1367.  URL https://doi.org/10.1214/15-AOS1367.
  * Zhong and Wang [2023] Qixian Zhong and Jane-Ling Wang.  Neural networks for partially linear quantile regression.  _Journal of Business & Economic Statistics_, 0(0):1–12, 2023.  doi: 10.1080/07350015.2023.2208183.  URL https://doi.org/10.1080/07350015.2023.2208183.
  * He and Shi [1996] Xuming He and Peide Shi.  Bivariate tensor-product b-splines in a partly linear model.  _Journal of Multivariate Analysis_ , 58(2):162–181, 1996.
  * Vansteelandt et al. [2022] Stijn Vansteelandt, Oliver Dukes, Kelly Van Lancker, and Torben Martinussen.  Assumption-Lean Cox Regression.  _Journal of the American Statistical Association_ , pages 1–10, 2022.  doi: 10.1080/01621459.2022.2126362.  URL https://doi.org/10.1080/01621459.2022.2126362.
  * van der Laan and Rose [2011] Mark J van der Laan and Sherri Rose.  _Targeted learning: causal inference for observational and experimental data_ , volume 4.  Springer Series in Statistics, New York, 2011.
  * Pfanzagl [1990] Johann Pfanzagl.  _Estimation in semiparametric models_ , pages 17–22.  Springer US, New York, NY, 1990.  ISBN 978-1-4612-3396-1.  doi: 10.1007/978-1-4612-3396-1\\_5.  URL https://doi.org/10.1007/978-1-4612-3396-1_5.
  * Bickel et al. [1993] Peter J Bickel, Chris AJ Klaassen, Peter J Bickel, Ya’acov Ritov, J Klaassen, Jon A Wellner, and YA’Acov Ritov.  _Efficient and adaptive estimation for semiparametric models_ , volume 4.  Baltimore: John Hopkins University Press, 1993.
  * Zheng and van der Laan [2011] Wenjing Zheng and Mark J. van der Laan.  _Cross-Validated Targeted Minimum-Loss-Based Estimation_ , pages 459–474.  Springer New York, New York, NY, 2011.  ISBN 978-1-4419-9782-1.  doi: 10.1007/978-1-4419-9782-1\\_27.  URL https://doi.org/10.1007/978-1-4419-9782-1_27.
  * van der Vaart [1998] A. W. van der Vaart.  _Functional Delta Method_ , page 291–303.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 1998.
  * van der Laan and Rubin [2006] Mark J. van der Laan and Daniel Rubin.  Targeted Maximum Likelihood Learning.  _The International Journal of Biostatistics_ , 2(1):Article 11, 2006.  doi: 10.2202/1557-4679.1043.  URL https://doi.org/10.2202/1557-4679.1043.
  * Yadlowsky [2022] Steve Yadlowsky.  Explaining practical differences between treatment effect estimators with high dimensional asymptotics.  _arXiv preprint arXiv:2203.12538_ , 2022.
  * Rosenbaum and Rubin [1983] Paul R. Rosenbaum and Donald B. Rubin.  The central role of the propensity score in observational studies for causal effects.  _Biometrika_ , 70(1):41–55, 04 1983.  ISSN 0006-3444.  doi: 10.1093/biomet/70.1.41.  URL https://doi.org/10.1093/biomet/70.1.41.
  * Rosenbaum [1989] Paul R. Rosenbaum.  Optimal matching for observational studies.  _Journal of the American Statistical Association_ , 84(408):1024–1032, 1989.  doi: 10.1080/01621459.1989.10478868.  URL https://www.tandfonline.com/doi/abs/10.1080/01621459.1989.10478868.
  * Sun et al. [2021] Shuo Sun, Erica E. M. Moodie, and Johanna G. Neslehová.  Causal inference for quantile treatment effects.  _Environmetrics_ , 32(4):e2668, 2021.  doi: https://doi.org/10.1002/env.2668.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/env.2668.
  * Athey et al. [2019] Susan Athey, Julie Tibshirani, and Stefan Wager.  Generalized random forests.  _The Annals of Statistics_ , 47(2):1148–1178, 2019.  doi: 10.1214/18-AOS1709.  URL https://doi.org/10.1214/18-AOS1709.
  * van der Laan et al. [2007] Mark J. van der Laan, Eric C Polley, and Alan E. Hubbard.  Super Learner.  _Statistical Applications in Genetics and Molecular Biology_ , 6(1):Article 25, 2007.  doi: doi:10.2202/1544-6115.1309.  URL https://doi.org/10.2202/1544-6115.1309.
  * Hofmeyr [2021] David P. Hofmeyr.  Fast Exact Evaluation of Univariate Kernel Sums.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 43(2):447–458, 2021.  doi: 10.1109/TPAMI.2019.2930501.
  * Hofmeyr [2022] David P. Hofmeyr.  Fast Kernel Smoothing in R with Applications to Projection Pursuit.  _Journal of Statistical Software_ , 101(3):1–33, 2022.  doi: 10.18637/jss.v101.i03.
  * Hejazi et al. [2022a] Nima S Hejazi, David Benkeser, Iván Díaz, and Mark J van der Laan.  Efficient estimation of modified treatment policy effects based on the generalized propensity score.  2022a.  URL https://arxiv.org/abs/2205.05777.
  * Hejazi et al. [2022b] Nima S Hejazi, David C Benkeser, and Mark J van der Laan.  haldensify: Highly adaptive lasso conditional density estimation.  https://github.com/nhejazi/haldensify, 2022b.  URL https://doi.org/10.5281/zenodo.3698329.  R package version 0.2.5.
  * Hejazi et al. [2022c] Nima S Hejazi, Mark J van der Laan, and David C Benkeser.  haldensify: Highly adaptive lasso conditional density estimation in R.  _Journal of Open Source Software_ , 2022c.  doi: 10.21105/joss.04522.  URL https://doi.org/10.21105/joss.04522.
  * Benkeser and Van Der Laan [2016] David Benkeser and Mark Van Der Laan.  The Highly Adaptive Lasso Estimator.  In _2016 IEEE International Conference on Data Science and Advanced Analytics (DSAA)_ , pages 689–696, 2016.  doi: 10.1109/DSAA.2016.93.
  * van der Laan [2017] Mark J van der Laan.  A Generally Efficient Targeted Minimum Loss Based Estimator based on the Highly Adaptive Lasso.  _The International Journal of Biostatistics_ , 13(2):20150097, 2017.  doi: doi:10.1515/ijb-2015-0097.  URL https://doi.org/10.1515/ijb-2015-0097.
  * Belloni et al. [2013a] Alexandre Belloni, Victor Chernozhukov, and Ying Wei.  Honest confidence regions for a regression parameter in logistic regression with a large number of controls.  cemmap working paper CWP67/13, Centre for Microdata Methods and Practice (cemmap), London, 2013a.  URL http://hdl.handle.net/10419/97416.
  * Gorasso et al. [2022] Vanessa Gorasso, Isabelle Moyersoen, Johan Van der Heyden, Karin De Ridder, Stefanie Vandevijvere, Stijn Vansteelandt, Delphine De Smedt, and Brecht Devleesschauwer.  Health care costs and lost productivity costs related to excess weight in Belgium.  _BMC Public Health_ , 22(1):1–11, 2022.
  * Belloni et al. [2013b] Alexandre Belloni, Victor Chernozhukov, and Christian Hansen.  Inference on Treatment Effects after Selection among High-Dimensional Controls†.  _The Review of Economic Studies_ , 81(2):608–650, 11 2013b.  ISSN 0034-6527.  doi: 10.1093/restud/rdt044.  URL https://doi.org/10.1093/restud/rdt044.
  * Robins [1986] James Robins.  A new approach to causal inference in mortality studies with a sustained exposure period—application to control of the healthy worker survivor effect.  _Mathematical Modelling_ , 7(9-12):1393–1512, 1986.
  * Demarest et al. [2013] Stefaan Demarest, Johan Van der Heyden, Rana Charafeddine, Sabine Drieskens, Lydia Gisle, and Jean Tafforeau.  Methodological basics and evolution of the Belgian health interview survey 1997–2008.  _Archives of Public Health_ , 71:1–10, 2013.
  * Huling et al. [2023] Jared D Huling, Noah Greifer, and Guanhuan Chen.  Independence weights for causal inference with continuous treatments.  _Journal of the American Statistical Association_ , 0(0):1–14, 2023.  doi: 10.1080/01621459.2023.2213485.  URL https://doi.org/10.1080/01621459.2023.2213485.
  * Chernozhukov et al. [2020] Victor Chernozhukov, Whitney Newey, Rahul Singh, and Vasilis Syrgkanis.  Adversarial estimation of Riesz representers.  _arXiv preprint arXiv:2101.00009_ , 2020.
  * Chernozhukov et al. [2022] Victor Chernozhukov, Whitney K Newey, and Rahul Singh.  Automatic debiased machine learning of causal and structural effects.  _Econometrica_ , 90(3):967–1027, 2022.
  * Guo and Shah [2023] F Richard Guo and Rajen D Shah.  Rank-transformed subsampling: inference for multiple data splitting and exchangeable p-values.  _arXiv preprint arXiv:2301.02739_ , 2023.
  * Hines et al. [2022] Oliver Hines, Oliver Dukes, Karla Diaz-Ordaz, and Stijn Vansteelandt.  Demystifying Statistical Learning Based on Efficient Influence Functions.  _The American Statistician_ , 76(3):292–304, 2022.  doi: 10.1080/00031305.2021.2021984.  URL https://doi.org/10.1080/00031305.2021.2021984.

Appendices

##  Appendix A Appendix A

###  A.1 Proof that (2) reduces to $\beta_{\tau}$ if model (1) holds.

Model (1) states that $Q_{\tau}(Y|A,L)=\beta_{\tau}A+\omega_{\tau}(L)$.
Substituting $Q_{\tau}(Y|A,L)$ and $Q_{\tau}(Y|A^{*},L)$ in (2) leads to

$\displaystyle\Psi_{\tau}$
$\displaystyle=\frac{E[(A-A^{*})\\{Q_{\tau}(Y|A,L)-Q_{\tau}(Y|A^{*},L)\\}]}{E\\{(A-A^{*})^{2}\\}}$
$\displaystyle=\frac{E[(A-A^{*})\\{\beta_{\tau}A+\omega_{\tau}(L)-\beta_{\tau}A^{*}+\omega_{\tau}(L)\\}]}{E\\{(A-A^{*})^{2}\\}}$
$\displaystyle=\frac{E[(A-A^{*})\\{\beta_{\tau}(A-A^{*})\\}]}{E\\{(A-A^{*})^{2}\\}}$
$\displaystyle=\beta_{\tau}.$

###  A.2 Proof that the estimands (2) and (3) are equivalent

We start from the numerator of (2):

$\displaystyle E[(A-A^{*})\\{Q_{\tau}(Y|A,L)-Q_{\tau}(Y|A^{*},L)\\}]$
$\displaystyle=E[\\{A-E(A|L)+E(A|L)-A^{*}\\}\\{Q_{\tau}(Y|A,L)-Q_{\tau}(Y|A^{*},L)\\}]$
$\displaystyle=E[\\{A-E(A|L)\\}\\{Q_{\tau}(Y|A,L)-Q_{\tau}(Y|A^{*},L)\\}]+E[\\{A^{*}-E(A|L)\\}\\{Q_{\tau}(Y|A^{*},L)-Q_{\tau}(Y|A,L)\\}]$
$\displaystyle=2E[\\{A-E(A|L)\\}Q_{\tau}(Y|A,L)]-E[\\{A-E(A|L)\\}Q_{\tau}(Y|A^{*},L)]-E[\\{A^{*}-E(A^{*}|L)\\}Q_{\tau}(Y|A,L)].$

We made use of the fact that $A$ and $A^{*}$ follow the same distribution
conditional on $L$. Moreover, the last two terms are zero because of the
(assumed) conditional independence of $A$ and $A^{*}$ given $L$:

$\displaystyle E[\\{A-E(A|L)\\}Q_{\tau}(Y|A^{*},L)]$
$\displaystyle=E\big{(}E[\\{A-E(A|L)\\}Q_{\tau}(Y|A^{*},L)|L]\big{)}$
$\displaystyle=E\big{(}E[\\{A-E(A|L)\\}|L]E\\{Q_{\tau}(Y|A^{*},L)|L\\}\big{)}$
$\displaystyle=E\big{[}0E\\{Q_{\tau}(Y|A^{*},L)|L\\}\big{]}$
$\displaystyle=0.$

The numerator of (3) reduces to

$E\left(\\{A-E(A|L)\\}\left[Q_{\tau}(Y|A,L)-E\\{Q_{\tau}(Y|A,L)|L\\}\right]\right)=E\left[\\{A-E(A|L)\\}Q_{\tau}(Y|A,L)\right].$

Thus, the numerator of (2) is equal to 2 times the numerator of (3). Now we
show that also the denominator of (2) is equal to 2 times the denominator of
(3), such that (2) and (3) are equivalent:

$\displaystyle E\\{(A-A^{*})^{2}\\}$
$\displaystyle=E[\\{A-E(A|L)+E(A|L)-A^{*}\\}^{2}]$
$\displaystyle=E[\\{A-E(A|L)\\}^{2}+\\{E(A|L)-A^{*}\\}^{2}-2\\{A-E(A|L)\\}\\{A^{*}-E(A|L)\\}]$
$\displaystyle=2E[\\{A-E(A|L)\\}^{2}]-2E[\\{A-E(A|L)\\}\\{A^{*}-E(A|L)\\}]$
$\displaystyle=2E[\\{A-E(A|L)\\}^{2}]-2E\big{(}E[\\{A-E(A|L)\\}\\{A^{*}-E(A|L)\\}|L]\big{)}$
$\displaystyle=2E[\\{A-E(A|L)\\}^{2}]-2E\big{(}E[\\{A-E(A|L)\\}|L]E[\\{A^{*}-E(A|L)\\}|L]\big{)}$
$\displaystyle=2E[\\{A-E(A|L)\\}^{2}].$

##  Appendix B Appendix B

###  B.1 Proof of Theorem 1

First note that estimand (3) can be equivalently written as

$\Psi_{\tau}=\frac{E[\\{A-E(A|L)\\}Q_{\tau}(Y|A,L)]}{E[\\{A-E(A|L)\\}^{2}]}.$
(25)

We will use (25) to calculate the efficient influence function of (3).
Furthermore, we follow the strategy described in Hines et al. [2022] as well
as the notation. Therefore, consider the one-dimensional parametric submodel
$\mathcal{P}_{t}=(1-t)\mathcal{P}+t\widetilde{\mathcal{P}}$, where
$\mathcal{P}$ is the true data-generating law and $\widetilde{\mathcal{P}}$ is
a point mass data generating distribution, such that the density of
$\mathcal{P}_{t}$ is given by
$f_{t}(y,a,l)=(1-t)f(y,a,l)+t\mathbbm{1}_{(\tilde{y},\tilde{a};\tilde{l})}(y,a,l)$,
where $\mathbbm{1}_{(\tilde{y},\tilde{a};\tilde{l})}(y,a,l)$ is a Dirac delta
function and $f(y,a,l)$ is the joint density of $Y$, $A$ and $L$. In the same
way, we also define all other functions, e.g.
$f_{t}(a|l)=f_{t}(a,l)/f_{t}(l)$. To calculate the efficient influence
function, we first give some intermediate results.

######  Lemma 1.

$\displaystyle\frac{d}{dt}Q_{\tau,t}(Y|A,L)\Big{\lvert}_{t=0}$
$\displaystyle=\frac{\mathbbm{1}_{\tilde{a},\tilde{l}}(A,L)}{f\\{Q_{\tau}(Y|A,L),A,L\\}}\Big{[}\tau-I\\{\tilde{y}\leq
Q_{\tau}(Y|A,L)\\}\Big{]}$

###### Proof.

Starting from

$\tau=\int_{-\infty}^{Q_{\tau}(Y|A,L)}f(Y|A,L)\,dY=\int_{-\infty}^{Q_{\tau}(Y|A,L)}\frac{f(Y,A,L)}{f(A,L)}\,dY$

and using the Leibniz integral rule, one can find

$\displaystyle 0$
$\displaystyle=\frac{d}{dt}\left(\int_{-\infty}^{Q_{\tau,t}(Y|A,L)}\frac{f_{t}(Y,A,L)}{f_{t}(A,L)}\,dY\right)\Bigg{\lvert}_{t=0}$
$\displaystyle=\frac{f(Q_{\tau}(Y|A,L),A,L)}{f(A,L)}\frac{d}{dt}Q_{\tau,t}(Y|A,L)\Big{\lvert}_{t=0}+\int_{-\infty}^{Q_{\tau}(Y|A,L)}\frac{\mathbbm{1}_{\tilde{y},\tilde{a},\tilde{l}}(Y,A,L)-f(Y,A,L)}{f(A,L)}\,dY$
$\displaystyle\qquad\qquad\qquad-\int_{-\infty}^{Q_{\tau}(Y|A,L)}\frac{[\mathbbm{1}_{\tilde{a},\tilde{l}}(A,L)-f(A,L)]f(Y,A,L)}{f(A,L)^{2}}\,dY$
$\displaystyle=\frac{f(Q_{\tau}(Y|A,L),A,L)}{f(A,L)}\frac{d}{dt}Q_{\tau,t}(Y|A,L)\Big{\lvert}_{t=0}+\int_{-\infty}^{Q_{\tau}(Y|A,L)}\frac{\mathbbm{1}_{\tilde{y},\tilde{a},\tilde{l}}(Y,A,L)}{f(A,L)}\,dY$
$\displaystyle\qquad\qquad\qquad-\int_{-\infty}^{Q_{\tau}(Y|A,L)}\frac{\mathbbm{1}_{\tilde{a},\tilde{l}}(A,L)f(Y,A,L)}{f(A,L)^{2}}\,dY$
$\displaystyle=\frac{f(Q_{\tau}(Y|A,L),A,L)}{f(A,L)}\frac{d}{dt}Q_{\tau,t}(Y|A,L)\Big{\lvert}_{t=0}+\frac{\mathbbm{1}_{\tilde{a},\tilde{l}}(A,L)}{f(A,L)}\int_{-\infty}^{Q_{\tau}(Y|A,L)}\mathbbm{1}_{\tilde{y}}(Y)\,dY$
$\displaystyle\qquad\qquad\qquad-\frac{\mathbbm{1}_{\tilde{a},\tilde{l}}(A,L)}{f(A,L)}\int_{-\infty}^{Q_{\tau}(Y|A,L)}f(Y|A,L)\,dY$
$\displaystyle=\frac{f(Q_{\tau}(Y|A,L),A,L)}{f(A,L)}\frac{d}{dt}Q_{\tau,t}(Y|A,L)\Big{\lvert}_{t=0}+\frac{\mathbbm{1}_{\tilde{a},\tilde{l}}(A,L)}{f(A,L)}I\\{\tilde{y}\leq
Q_{\tau}(Y|A,L)\\}$
$\displaystyle\qquad\qquad\qquad-\frac{\mathbbm{1}_{\tilde{a},\tilde{l}}(A,L)}{f(A,L)}\tau.$

This leads to

$\frac{d}{dt}Q_{\tau,t}(Y|A,L)\Big{\lvert}_{t=0}=\frac{\mathbbm{1}_{\tilde{a},\tilde{l}}(A,L)}{f(Q_{\tau}(Y|A,L),A,L)}\Big{(}\tau-I\\{\tilde{y}\leq
Q_{\tau}(Y|A,L)\\}\Big{)}.$

∎

######  Lemma 2.

$\displaystyle\frac{d}{dt}(A-E_{t}(A|L))Q_{\tau,t}(Y|A,L)\Big{\lvert}_{t=0}$
$\displaystyle=-\frac{\mathbbm{1}_{\tilde{l}}(L)}{f(L)}(\tilde{a}-E(A|L))Q_{\tau}(Y|A,L)$
$\displaystyle\qquad\qquad+\\{A-E(A|L)\\}\mathbbm{1}_{\tilde{a},\tilde{l}}(A,L)\frac{\tau-I\\{\tilde{y}\leq
Q_{\tau}(Y|A,L)\\}}{f(Q_{\tau}(Y|A,L),A,L)}$

###### Proof.

$\displaystyle\frac{d}{dt}(A-E_{t}(A|L))Q_{\tau,t}(Y|A,L)\Big{\lvert}_{t=0}$
$\displaystyle=\frac{d}{dt}(A-E_{t}(A|L))\Big{\lvert}_{t=0}Q_{\tau}(Y|A,L)$
$\displaystyle\qquad\qquad+(A-E(A|L))\frac{d}{dt}Q_{\tau,t}(Y|A,L)\Big{\lvert}_{t=0}$
$\displaystyle=-\frac{\mathbbm{1}_{\tilde{l}}(L)}{f(L)}(\tilde{a}-E(A|L))Q_{\tau}(Y|A,L)$
$\displaystyle\qquad\qquad+\\{A-E(A|L)\\}\mathbbm{1}_{\tilde{a},\tilde{l}}(A,L)\frac{\tau-I\\{\tilde{y}\leq
Q_{\tau}(Y|A,L)\\}}{f(Q_{\tau}(Y|A,L),A,L)}$

∎

######  Lemma 3.

$\displaystyle\frac{d}{dt}(A-E_{t}(A|L))^{2}\Big{\lvert}_{t=0}$
$\displaystyle=-2\\{A-E(A|L)\\}(\tilde{a}-E(A|L))\frac{\mathbbm{1}_{\tilde{l}}(L)}{f(L)}$

###### Proof.

$\displaystyle\frac{d}{dt}(A-E_{t}(A|L))^{2}\Big{\lvert}_{t=0}$
$\displaystyle=-2\\{A-E(A|L)\\}\frac{d}{dt}E_{t}(A|L)\Big{\lvert}_{t=0}$
$\displaystyle=-2\\{A-E(A|L)\\}(\tilde{a}-E(A|L))\frac{\mathbbm{1}_{\tilde{l}}(L)}{f(L)}$

∎

Taking the derivative of $\Psi_{\tau}(\mathcal{P}_{t})$ with respect to $t$,
evaluated at $t=0$, results in

$\displaystyle\frac{d}{dt}\Psi_{\tau}(\mathcal{P}_{t})\Big{\lvert}_{t=0}$
$\displaystyle=\frac{\frac{d}{dt}E_{t}[(A-E_{t}(A|L))Q_{\tau,t}(Y|A,L)]\Big{\lvert}_{t=0}}{E[\\{A-E(A|L)\\}^{2}]}$
$\displaystyle\qquad\qquad\qquad\qquad-\frac{\Psi_{\tau}}{E[\\{A-E(A|L)\\}^{2}]}\frac{d}{dt}E_{t}[(A-E_{t}(A|L))^{2}]\Big{\lvert}_{t=0}$
$\displaystyle=\frac{(\tilde{a}-E(A|L=\tilde{l}))Q_{\tau}(Y|A=\tilde{a},L=\tilde{l})}{E[\\{A-E(A|L)\\}^{2}]}-\Psi_{\tau}+\frac{E\left[\frac{d}{dt}(A-E_{t}(A|L))Q_{\tau,t}(Y|A,L)\Big{\lvert}_{t=0}\right]}{E[\\{A-E(A|L)\\}^{2}]}$
$\displaystyle\qquad-\frac{\Psi_{\tau}}{E[\\{A-E(A|L)\\}^{2}]}\Bigg{[}(\tilde{a}-E(A|L=\tilde{l}))^{2}-E[\\{A-E(A|L)\\}^{2}]+E\left(\frac{d}{dt}(A-E_{t}(A|L))^{2}\Big{\lvert}_{t=0}\right)\Bigg{]}$
$\displaystyle=\frac{1}{E[\\{A-E(A|L)\\}^{2}]}\Bigg{[}(\tilde{a}-E(A|L=\tilde{l}))Q_{\tau}(Y|A=\tilde{a},L=\tilde{l})$
$\displaystyle\qquad\qquad-\iint\frac{\mathbbm{1}_{\tilde{l}}(L)}{f(L)}(\tilde{a}-E(A|L))Q_{\tau}(Y|A,L)f(A,L)\,dA\,dL$
$\displaystyle\qquad\qquad+\iint\\{A-E(A|L)\\}\mathbbm{1}_{\tilde{a},\tilde{l}}(A,L)\frac{\tau-I\\{\tilde{y}\leq
Q_{\tau}(Y|A,L)\\}}{f(Q_{\tau}(Y|A,L),A,L)}f(A,L)\,dA\,dL$
$\displaystyle\qquad\qquad-\Psi_{\tau}\Bigg{(}(\tilde{a}-E(A|L=\tilde{l}))^{2}-2\iint\\{A-E(A|L)\\}(\tilde{a}-E(A|L))\frac{\mathbbm{1}_{\tilde{l}}(L)}{f(L)}f(A,L)\,dA\,dL\Bigg{)}\Bigg{]}$
$\displaystyle=\frac{1}{E[\\{A-E(A|L)\\}^{2}]}\Bigg{[}(\tilde{a}-E(A|L=\tilde{l}))Q_{\tau}(Y|A=\tilde{a},L=\tilde{l})$
$\displaystyle\qquad\qquad-(\tilde{a}-E(A|L=\tilde{l}))\int
Q_{\tau}(Y|A,L=\tilde{l})f(A|L=\tilde{l})\,dA$
$\displaystyle\qquad\qquad+(\tilde{a}-E(A|L=\tilde{l}))\frac{\tau-I\\{\tilde{y}\leq
Q_{\tau}(Y|A=\tilde{a},L=\tilde{l})\\}}{f(Q_{\tau}(Y|A=\tilde{a},L=\tilde{l})|A=\tilde{a},L=\tilde{l})}$
$\displaystyle\qquad\qquad-\Psi_{\tau}\Bigg{(}(\tilde{a}-E(A|L=\tilde{l}))^{2}-2(\tilde{a}-E(A|L=\tilde{l}))\int(A-E(A|L=\tilde{l}))f(A|L=\tilde{l})\,dA\Bigg{)}\Bigg{]}$
$\displaystyle=\frac{\tilde{a}-E(A|L=\tilde{l})}{E[\\{A-E(A|L)\\}^{2}]}\Bigg{[}Q_{\tau}(Y|A=\tilde{a},L=\tilde{l})-E\Big{(}Q_{\tau}(Y|A,L=\tilde{l})\big{\lvert}L=\tilde{l}\Big{)}$
$\displaystyle\qquad\qquad\qquad\qquad+\frac{\tau-I\\{\tilde{y}\leq
Q_{\tau}(Y|A=\tilde{a},L=\tilde{l})\\}}{f(Q_{\tau}(Y|A=\tilde{a},L=\tilde{l})|A=\tilde{a},L=\tilde{l})}$
$\displaystyle\qquad\qquad\qquad\qquad-\Psi_{\tau}\Bigg{(}(\tilde{a}-E(A|L=\tilde{l}))-2E\big{(}A-E(A|L=\tilde{l})|L=\tilde{l}\big{)}\Bigg{)}\Bigg{]}$
$\displaystyle=\frac{\tilde{a}-E(A|L=\tilde{l})}{E[\\{A-E(A|L)\\}^{2}]}\Bigg{[}Q_{\tau}(Y|A=\tilde{a},L=\tilde{l})-E\Big{(}Q_{\tau}(Y|A,L=\tilde{l})\big{\lvert}L=\tilde{l}\Big{)}$
$\displaystyle\qquad\qquad\qquad\qquad+\frac{\tau-I\\{\tilde{y}\leq
Q_{\tau}(Y|A=\tilde{a},L=\tilde{l})\\}}{f(Q_{\tau}(Y|A=\tilde{a},L=\tilde{l})|A=\tilde{a},L=\tilde{l})}-\Psi_{\tau}(\tilde{a}-E(A|L=\tilde{l}))\Bigg{]},$

where we use Lemma 2 and Lemma 3 in the third equation. Going back to the
random variables, this results in the following efficient influence function:

$\displaystyle\frac{A-E(A|L)}{E[\\{A-E(A|L)\\}^{2}]}\Bigg{[}Q_{\tau}(Y|A,L)-E\left\\{Q_{\tau}(Y|A,L)|L\right\\}+\frac{\tau-I\\{Y\leq
Q_{\tau}(Y|A,L)\\}}{f\\{Q_{\tau}(Y|A,L)|A,L\\}}-\Psi_{\tau}\\{A-E(A|L)\\}\Bigg{]}.$

###  B.2 Proof of Theorem 2

Using a von Mises expansion (see e.g. van der Vaart [1998]), we have

$\displaystyle\hat{\Psi}_{\tau}-\Psi_{\tau}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\frac{A_{i}-E(A_{i}|L_{i})}{E[\\{A-E(A|L)\\}^{2}]}\Bigg{[}Q_{\tau}(Y_{i}|A_{i},L_{i})-E\left\\{Q_{\tau}(Y|A,L)|L\right\\}$
$\displaystyle\qquad\qquad+\frac{\tau-I\\{Y_{i}\leq
Q_{\tau}(Y_{i}|A_{i},L_{i})\\}}{f\\{Q_{\tau}(Y_{i}|A_{i},L_{i})|A_{i},L_{i}\\}}-\Psi_{\tau}(A_{i}-E(A_{i}|L_{i}))\Bigg{]}+R_{1}+R_{2}.$

Here, the first term is $O_{p}(n^{-1/2})$ and

$\displaystyle R_{1}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\Bigg{[}\frac{A_{i}-\hat{E}(A_{i}|L_{i})}{\hat{E}[(A_{i}-\hat{E}(A_{i}|L_{i}))^{2}]}\Bigg{\\{}\hat{Q}_{\tau}(Y_{i}|A_{i},L_{i})-\hat{E}\Big{(}\hat{Q}_{\tau}(Y_{i}|A_{i},L_{i})\big{\lvert}L_{i}\Big{)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad+\frac{\tau-I\\{Y_{i}\leq\hat{Q}_{\tau}(Y_{i}|A_{i},L_{i})\\}}{\hat{f}\\{\hat{Q}_{\tau}(Y_{i}|A_{i},L_{i})|A_{i},L_{i}\\}}-\hat{\Psi}_{\tau}(A_{i}-\hat{E}(A_{i}|L_{i}))\Bigg{\\}}$
$\displaystyle\qquad\qquad\qquad-\frac{A_{i}-E(A_{i}|L_{i})}{E[\\{A-E(A|L)\\}^{2}]}\Bigg{\\{}Q_{\tau}(Y_{i}|A_{i},L_{i})-E\left\\{Q_{\tau}(Y|A,L)|L\right\\}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad+\frac{\tau-I\\{Y_{i}\leq
Q_{\tau}(Y_{i}|A_{i},L_{i})\\}}{f\\{Q_{\tau}(Y_{i}|A_{i},L_{i})|A_{i},L_{i}\\}}-\Psi_{\tau}(A_{i}-E(A_{i}|L_{i}))\Bigg{\\}}\Bigg{]}$
$\displaystyle\qquad-E\Bigg{[}\frac{A-\hat{E}(A|L)}{\hat{E}[(A-\hat{E}(A|L))^{2}]}\Bigg{\\{}\hat{Q}_{\tau}(Y|A,L)-\hat{E}\Big{(}\hat{Q}_{\tau}(Y|A,L)\big{\lvert}L\Big{)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad+\frac{\tau-I\\{Y\leq\hat{Q}_{\tau}(Y|A,L)\\}}{\hat{f}(\hat{Q}_{\tau}(Y|A,L)|A,L)}-\hat{\Psi}_{\tau}(A-\hat{E}(A|L))\Bigg{\\}}$
$\displaystyle\qquad\qquad\qquad-\frac{A-E(A|L)}{E[\\{A-E(A|L)\\}^{2}]}\Bigg{\\{}Q_{\tau}(Y|A,L)-E\left\\{Q_{\tau}(Y|A,L)|L\right\\}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad+\frac{\tau-I\\{Y\leq
Q_{\tau}(Y|A,L)\\}}{f\\{Q_{\tau}(Y|A,L)|A,L\\}}-\Psi_{\tau}\\{A-E(A|L)\\}\Bigg{\\}}\Bigg{]}$

and

$\displaystyle R_{2}$
$\displaystyle=\hat{\Psi}_{\tau}-\Psi_{\tau}+E\Bigg{[}\frac{A-\hat{E}(A|L)}{\hat{E}[(A-\hat{E}(A|L))^{2}]}\Bigg{\\{}\hat{Q}_{\tau}(Y|A,L)-\hat{E}\Big{(}\hat{Q}_{\tau}(Y|A,L)\big{\lvert}L\Big{)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad+\frac{\tau-I\\{Y\leq\hat{Q}_{\tau}(Y|A,L)\\}}{\hat{f}(\hat{Q}_{\tau}(Y|A,L)|A,L)}-\hat{\Psi}_{\tau}(A-\hat{E}(A|L))\Bigg{\\}}\Bigg{]}.$

The term $R_{1}$ can be shown to be $o_{p}(n^{-1/2})$ by using sample-
splitting. A detailed elaboration of this can be found along the lines of the
Supplementary Material of Vansteelandt and Dukes [2022]. In what follows we
take a deeper look into $R_{2}$ to understand under what conditions it is
$o_{p}(n^{-1/2})$. From the definition of $R_{2}$, it follows that

$\displaystyle R_{2}$
$\displaystyle=\hat{\Psi}_{\tau}-\Psi_{\tau}-\hat{\Psi}_{\tau}\frac{E[(A-\hat{E}(A|L))^{2}]}{\hat{E}[(A-\hat{E}(A|L))^{2}]}$
$\displaystyle\qquad\qquad+E\Bigg{[}\frac{A-\hat{E}(A|L)}{\hat{E}[(A-\hat{E}(A|L))^{2}]}\Bigg{\\{}\hat{Q}_{\tau}(Y|A,L)-\hat{E}\Big{(}\hat{Q}_{\tau}(Y|A,L)\big{\lvert}L\Big{)}+\frac{\tau-I\\{Y\leq\hat{Q}_{\tau}(Y|A,L)\\}}{\hat{f}(\hat{Q}_{\tau}(Y|A,L)|A,L)}\Bigg{\\}}\Bigg{]}$
$\displaystyle=(\hat{\Psi}_{\tau}-\Psi_{\tau})\Bigg{[}1-\frac{E[(A-\hat{E}(A|L))^{2}]}{\hat{E}[(A-\hat{E}(A|L))^{2}]}\Bigg{]}$
$\displaystyle\qquad\qquad+E\Bigg{[}\frac{A-\hat{E}(A|L)}{\hat{E}[(A-\hat{E}(A|L))^{2}]}\Bigg{\\{}\hat{Q}_{\tau}(Y|A,L)-\hat{E}\Big{(}\hat{Q}_{\tau}(Y|A,L)\big{\lvert}L\Big{)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad+\frac{\tau-I\\{Y\leq\hat{Q}_{\tau}(Y|A,L)\\}}{\hat{f}(\hat{Q}_{\tau}(Y|A,L)|A,L)}-\Psi_{\tau}(A-\hat{E}(A|L))\Bigg{\\}}\Bigg{]}.$

Here, the first term is $o_{p}(|\hat{\Psi}_{\tau}-\Psi_{\tau}|)$, and thus a
lower order term. To understand the behaviour of the remaining terms, we make
use of

$\displaystyle\Psi_{\tau}$
$\displaystyle=E\Bigg{[}\frac{A-E(A|L)}{E[(A-E(A|L))^{2}]}\Bigg{\\{}Q_{\tau}(Y|A,L)-E\Big{(}Q_{\tau}(Y|A,L)\big{\lvert}L\Big{)}+\frac{\tau-I\\{Y\leq
Q_{\tau}(Y|A,L)\\}}{f\\{Q_{\tau}(Y|A,L)|A,L\\}}\Bigg{\\}}\Bigg{]}.$

To simplify the notation in the following calculations, let
$\zeta=Q_{\tau}(Y|A,L)$. Then the remaining terms can be written as

$\displaystyle\frac{1}{\hat{E}[(A-\hat{E}(A|L))^{2}]}\Bigg{(}E\left[(A-\hat{E}(A|L))\Big{\\{}\hat{\zeta}-\hat{E}(\hat{\zeta}|L)+\frac{\tau-I\\{Y\leq\hat{\zeta}\\}}{\hat{f}(\hat{\zeta}|A,L)}\Big{\\}}\right]$
$\displaystyle\qquad\qquad-E\left[(A-E(A|L))\Big{\\{}\zeta-E(\zeta|L)+\frac{\tau-I\\{Y\leq\zeta\\}}{f(\zeta|A,L)}\Big{\\}}\right]\frac{E[(A-\hat{E}(A|L))^{2}]}{E[(A-E(A|L))^{2}]}\Bigg{)}$
$\displaystyle=\frac{1}{\hat{E}[(A-\hat{E}(A|L))^{2}]}\Bigg{(}E\left[(A-\hat{E}(A|L))\Big{\\{}\hat{\zeta}-\hat{E}(\hat{\zeta}|L)+\frac{\tau-P(Y\leq\hat{\zeta}|A,L)}{\hat{f}(\hat{\zeta}|A,L)}\Big{\\}}\right]$
$\displaystyle\qquad\qquad-E\left[(A-E(A|L))\Big{\\{}\zeta-E(\zeta|L)+\frac{\tau-P(Y\leq\zeta|A,L)}{f(\zeta|A,L)}\Big{\\}}\right]\frac{E[(A-\hat{E}(A|L))^{2}]}{E[(A-E(A|L))^{2}]}\Bigg{)}$
$\displaystyle=\frac{1}{\hat{E}[(A-\hat{E}(A|L))^{2}]}\Bigg{(}E\left[(A-\hat{E}(A|L))\Big{\\{}\hat{\zeta}-\hat{E}(\hat{\zeta}|L)+\frac{\tau-P(Y\leq\hat{\zeta}|A,L)}{\hat{f}(\hat{\zeta}|A,L)}\Big{\\}}\right]$
$\displaystyle\qquad\qquad-E\left[(A-E(A|L))\Big{\\{}\zeta-E(\zeta|L)\Big{\\}}\right]\frac{E[(A-\hat{E}(A|L))^{2}]}{E[(A-E(A|L))^{2}]}\Bigg{)}.$

Next, using a Taylor expansion around $\tau$, one can write

$P(Y\leq\hat{\zeta}|A,L)=\tau+f(\zeta|A,L)(\hat{\zeta}-\zeta)+O_{p}\\{(\hat{\zeta}-\zeta)^{2}\\}.$

This reduces the remaining terms to

$\displaystyle\frac{1}{\hat{E}[(A-\hat{E}(A|L))^{2}]}\Bigg{(}E\left[(A-\hat{E}(A|L))\Big{\\{}\hat{\zeta}-\hat{E}(\hat{\zeta}|L)-\frac{f(\zeta|A,L)}{\hat{f}(\hat{\zeta}|A,L)}(\hat{\zeta}-\zeta)+O_{p}\\{(\hat{\zeta}-\zeta)^{2}\\}\Big{\\}}\right]$
$\displaystyle\qquad\qquad-E\left[(A-E(A|L))\Big{\\{}\zeta-E(\zeta|L)\Big{\\}}\right]\frac{E[(A-\hat{E}(A|L))^{2}]}{E[(A-E(A|L))^{2}]}\Bigg{)}$
$\displaystyle=\frac{1}{\hat{E}[(A-\hat{E}(A|L))^{2}]}\Bigg{(}E\left[(A-\hat{E}(A|L))\Big{\\{}\zeta-\hat{E}(\hat{\zeta}|L)-\left(\frac{f(\zeta|A,L)}{\hat{f}(\hat{\zeta}|A,L)}-1\right)(\hat{\zeta}-\zeta)+O_{p}\\{(\hat{\zeta}-\zeta)^{2}\\}\Big{\\}}\right]$
$\displaystyle\qquad\qquad-E\left[(A-E(A|L))\Big{\\{}\zeta-E(\zeta|L)\Big{\\}}\right]\frac{E[(A-\hat{E}(A|L))^{2}]}{E[(A-E(A|L))^{2}]}\Bigg{)}$
$\displaystyle=\frac{1}{\hat{E}[(A-\hat{E}(A|L))^{2}]}\Bigg{(}E\left[(A-E(A|L))\Big{\\{}\zeta-\hat{E}(\hat{\zeta}|L)-\left(\frac{f(\zeta|A,L)}{\hat{f}(\hat{\zeta}|A,L)}-1\right)(\hat{\zeta}-\zeta)+O_{p}\\{(\hat{\zeta}-\zeta)^{2}\\}\Big{\\}}\right]$
$\displaystyle\qquad\qquad+E\left[(E(A|L)-\hat{E}(A|L))\Big{\\{}\zeta-\hat{E}(\hat{\zeta}|L)-\left(\frac{f(\zeta|A,L)}{\hat{f}(\hat{\zeta}|A,L)}-1\right)(\hat{\zeta}-\zeta)+O_{p}\\{(\hat{\zeta}-\zeta)^{2}\\}\Big{\\}}\right]$
$\displaystyle\qquad\qquad-E\left[(A-E(A|L))\Big{\\{}\zeta-E(\zeta|L)\Big{\\}}\right]$
$\displaystyle\qquad\qquad+\left(1-\frac{E[(A-\hat{E}(A|L))^{2}]}{E[(A-E(A|L))^{2}]}\right)E\left[(A-E(A|L))\Big{\\{}\zeta-E(\zeta|L)\Big{\\}}\right]\Bigg{)}$
$\displaystyle=\frac{1}{\hat{E}[(A-\hat{E}(A|L))^{2}]}\Bigg{(}E\left[(A-E(A|L))\Big{\\{}-\hat{E}(\hat{\zeta}|L)-\left(\frac{f(\zeta|A,L)}{\hat{f}(\hat{\zeta}|A,L)}-1\right)(\hat{\zeta}-\zeta)+O_{p}\\{(\hat{\zeta}-\zeta)^{2}\\}\Big{\\}}\right]$
$\displaystyle\qquad\qquad+E\left[(E(A|L)-\hat{E}(A|L))\Big{\\{}\zeta-\hat{E}(\hat{\zeta}|L)-\left(\frac{f(\zeta|A,L)}{\hat{f}(\hat{\zeta}|A,L)}-1\right)(\hat{\zeta}-\zeta)+O_{p}\\{(\hat{\zeta}-\zeta)^{2}\\}\Big{\\}}\right]$
$\displaystyle\qquad\qquad-E\left[(A-E(A|L))\Big{\\{}-E(\zeta|L)\Big{\\}}\right]$
$\displaystyle\qquad\qquad+\left(1-\frac{E[(A-\hat{E}(A|L))^{2}]}{E[(A-E(A|L))^{2}]}\right)E\left[(A-E(A|L))\Big{\\{}\zeta-E(\zeta|L)\Big{\\}}\right]\Bigg{)}$
$\displaystyle=\frac{1}{\hat{E}[(A-\hat{E}(A|L))^{2}]}\Bigg{(}E\left[(A-E(A|L))\Big{\\{}O_{p}\\{(\hat{\zeta}-\zeta)^{2}\\}-\left(\frac{f(\zeta|A,L)}{\hat{f}(\hat{\zeta}|A,L)}-1\right)(\hat{\zeta}-\zeta)\Big{\\}}\right]$
$\displaystyle\qquad+E\left[(E(A|L)-\hat{E}(A|L))\Big{\\{}\zeta-\hat{E}(\hat{\zeta}|L)+E(\zeta|L)-E(\zeta|L)-\left(\frac{f(\zeta|A,L)}{\hat{f}(\hat{\zeta}|A,L)}-1\right)(\hat{\zeta}-\zeta)+O_{p}\\{(\hat{\zeta}-\zeta)^{2}\\}\Big{\\}}\right]$
$\displaystyle\qquad\qquad+\left(1-\frac{E[(A-\hat{E}(A|L))^{2}]}{E[(A-E(A|L))^{2}]}\right)E\left[(A-E(A|L))\Big{\\{}\zeta-E(\zeta|L)\Big{\\}}\right]\Bigg{)}$
$\displaystyle=\frac{1}{\hat{E}[(A-\hat{E}(A|L))^{2}]}\Bigg{(}E\left[(A-E(A|L))\Big{\\{}O_{p}\\{(\hat{\zeta}-\zeta)^{2}\\}-\left(\frac{f(\zeta|A,L)}{\hat{f}(\hat{\zeta}|A,L)}-1\right)(\hat{\zeta}-\zeta)\Big{\\}}\right]$
$\displaystyle\qquad\qquad+E\left[(E(A|L)-\hat{E}(A|L))\Big{\\{}\zeta-E(\zeta|L)-\left(\frac{f(\zeta|A,L)}{\hat{f}(\hat{\zeta}|A,L)}-1\right)(\hat{\zeta}-\zeta)+O_{p}\\{(\hat{\zeta}-\zeta)^{2}\\}\Big{\\}}\right]$
$\displaystyle\qquad\qquad+E\left[(E(A|L)-\hat{E}(A|L))\Big{\\{}E(\zeta|L)-\hat{E}(\hat{\zeta}|L)\Big{\\}}\right]$
$\displaystyle\qquad\qquad+\left(1-\frac{E[(A-\hat{E}(A|L))^{2}]}{E[(A-E(A|L))^{2}]}\right)E\left[(A-E(A|L))\Big{\\{}\zeta-E(\zeta|L)\Big{\\}}\right]\Bigg{)}$
$\displaystyle=\frac{1}{\hat{E}[(A-\hat{E}(A|L))^{2}]}\Bigg{(}E\left[(A-E(A|L))\Big{\\{}O_{p}\\{(\hat{\zeta}-\zeta)^{2}\\}-\left(\frac{f(\zeta|A,L)}{\hat{f}(\hat{\zeta}|A,L)}-1\right)(\hat{\zeta}-\zeta)\Big{\\}}\right]$
$\displaystyle\qquad\qquad+E\left[(E(A|L)-\hat{E}(A|L))\Big{\\{}O_{p}\\{(\hat{\zeta}-\zeta)^{2}\\}-\left(\frac{f(\zeta|A,L)}{\hat{f}(\hat{\zeta}|A,L)}-1\right)(\hat{\zeta}-\zeta)\Big{\\}}\right]$
$\displaystyle\qquad\qquad+E\left[(E(A|L)-\hat{E}(A|L))\Big{\\{}E(\zeta|L)-\hat{E}(\hat{\zeta}|L)\Big{\\}}\right]$
$\displaystyle\qquad\qquad+\left(1-\frac{E[(A-E(A|L)+E(A|L)-\hat{E}(A|L))^{2}]}{E[(A-E(A|L))^{2}]}\right)E\left[(A-E(A|L))\Big{\\{}\zeta-E(\zeta|L)\Big{\\}}\right]\Bigg{)}$
$\displaystyle=\frac{1}{\hat{E}[(A-\hat{E}(A|L))^{2}]}\Bigg{(}E\left[(A-\hat{E}(A|L))\Big{\\{}O_{p}\\{(\hat{\zeta}-\zeta)^{2}\\}-\left(\frac{f(\zeta|A,L)}{\hat{f}(\hat{\zeta}|A,L)}-1\right)(\hat{\zeta}-\zeta)\Big{\\}}\right]$
$\displaystyle\qquad\qquad+E\left[(E(A|L)-\hat{E}(A|L))\Big{\\{}E(\zeta|L)-\hat{E}(\hat{\zeta}|L)\Big{\\}}\right]$
$\displaystyle\qquad\qquad-\frac{E[2(A-E(A|L))(E(A|L)-\hat{E}(A|L))+(E(A|L)-\hat{E}(A|L))^{2}]}{E[(A-E(A|L))^{2}]}E\left[(A-E(A|L))\Big{\\{}\zeta-E(\zeta|L)\Big{\\}}\right]\Bigg{)}$
$\displaystyle=\frac{1}{\hat{E}[(A-\hat{E}(A|L))^{2}]}\Bigg{(}E\left[(A-\hat{E}(A|L))O_{p}\\{(\hat{\zeta}-\zeta)^{2}\\}\right]+E\left[(A-\hat{E}(A|L))\left(1-\frac{f(\zeta|A,L)}{\hat{f}(\hat{\zeta}|A,L)}\right)(\hat{\zeta}-\zeta)\right]$
$\displaystyle\qquad\qquad+E\left[(E(A|L)-\hat{E}(A|L))\Big{\\{}E(\zeta|L)-\hat{E}(\hat{\zeta}|L)\Big{\\}}\right]$
$\displaystyle\qquad\qquad-\frac{E[(E(A|L)-\hat{E}(A|L))^{2}]}{E[(A-E(A|L))^{2}]}E\left[(A-E(A|L))\Big{\\{}\zeta-E(\zeta|L)\Big{\\}}\right]\Bigg{)}$
$\displaystyle=\frac{1}{\hat{E}[(A-\hat{E}(A|L))^{2}]}\Bigg{(}E\left[(A-\hat{E}(A|L))O_{p}\\{(\hat{\zeta}-\zeta)^{2}\\}\right]+E\left[(A-\hat{E}(A|L))\left(1-\frac{f(\zeta|A,L)}{\hat{f}(\hat{\zeta}|A,L)}\right)(\hat{\zeta}-\zeta)\right]$
$\displaystyle\qquad\qquad+E\left[(E(A|L)-\hat{E}(A|L))\Big{\\{}E(\zeta|L)-\hat{E}(\hat{\zeta}|L)\Big{\\}}\right]-E[(E(A|L)-\hat{E}(A|L))^{2}]\Psi_{\tau}\Bigg{)}.$

Assuming $A-\hat{E}(A|L)=O_{p}(1)$ and using the Cauchy-Schwarz inequality,
the term $R_{2}$ is $o_{p}(n^{-1/2})$ under the following conditions:

$\displaystyle E\left[(\hat{Q}_{\tau}(Y|A,L)-Q_{\tau}(Y|A,L))^{2}\right]$
$\displaystyle=o_{p}(n^{-1/2}),$ $\displaystyle
E\left[\left(1-\frac{f\\{Q_{\tau}(Y|A,L)|A,L\\}}{\hat{f}(\hat{Q}_{\tau}(Y|A,L)|A,L)}\right)^{2}\right]^{1/2}E\left[(\hat{Q}_{\tau}(Y|A,L)-Q_{\tau}(Y|A,L))^{2}\right]^{1/2}$
$\displaystyle=o_{p}(n^{-1/2}),$ $\displaystyle
E\left[(E(A|L)-\hat{E}(A|L))^{2}\right]^{1/2}E\left[\left(E(Q_{\tau}(Y|A,L)|L)-\hat{E}(\hat{Q}_{\tau}(Y|A,L)|L)\right)^{2}\right]^{1/2}$
$\displaystyle=o_{p}(n^{-1/2}),$ $\displaystyle
E\left[(E(A|L)-\hat{E}(A|L))^{2}\right]$ $\displaystyle=o_{p}(n^{-1/2}).$

Note that if $\Psi_{\tau}=0$, the final condition is redundant.

###  B.3 Sketch of a proof that the sum in (11) is $o_{P}(n^{-1/2})$

Here, we give a sketch of a proof that

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left(A_{i}-\hat{E}(A_{i}|L_{i})\right)\Bigg{[}\frac{\tau-I\\{Y_{i}\leq\widetilde{Q}_{\tau}(Y_{i}|A_{i},L_{i})\\}}{\hat{f}(\widetilde{Q}_{\tau}(Y_{i}|A_{i},L_{i})|A_{i},L_{i})}\Bigg{]}=o_{P}(n^{-1/2}).$

We consider

$\widetilde{Q}_{\tau}(Y|A,L)=\hat{Q}_{\tau}(Y|A,L)+\epsilon
w(\hat{Q}_{\tau,i}),$

with $w(\hat{Q}_{\tau,i})$ defined as in (12). Thus, the aim is to prove that

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}w(\hat{Q}_{\tau,i})\left[\tau-I\\{Y_{i}\leq\hat{Q}_{\tau}(Y_{i}|A_{i},L_{i})+\epsilon
w(\hat{Q}_{\tau,i})\\}\right]$ (26)

is $o_{P}(n^{-1/2})$.

One can control the sign of (26) by changing the value of $\epsilon$. Indeed,
for $\epsilon\to+\infty$ and $w(\hat{Q}_{\tau,i})\geq 0$, we have

$w(\hat{Q}_{\tau,i})\left[\tau-I\\{Y_{i}\leq\hat{Q}_{\tau}(Y_{i}|A_{i},L_{i})+\epsilon
w(\hat{Q}_{\tau,i})\\}\right]=w(\hat{Q}_{\tau,i})(\tau-1)\leq 0.$

For $\epsilon\to+\infty$ and $w(\hat{Q}_{\tau,i})\leq 0$, we have

$w(\hat{Q}_{\tau,i})\left[\tau-I\\{Y_{i}\leq\hat{Q}_{\tau}(Y_{i}|A_{i},L_{i})+\epsilon
w(\hat{Q}_{\tau,i})\\}\right]=w(\hat{Q}_{\tau,i})\tau\leq 0,$

such that

$\lim_{\epsilon\to+\infty}\frac{1}{n}\sum_{i=1}^{n}w(\hat{Q}_{\tau,i})\left[\tau-I\\{Y_{i}\leq\hat{Q}_{\tau}(Y_{i}|A_{i},L_{i})+\epsilon
w(\hat{Q}_{\tau,i})\\}\right]\leq 0.$

Analogously, one can show

$\lim_{\epsilon\to-\infty}\frac{1}{n}\sum_{i=1}^{n}w(\hat{Q}_{\tau,i})\left[\tau-I\\{Y_{i}\leq\hat{Q}_{\tau}(Y_{i}|A_{i},L_{i})+\epsilon
w(\hat{Q}_{\tau,i})\\}\right]\geq 0.$

By slightly adjusting the value of $\epsilon$, terms from the summation can
change from $w(\hat{Q}_{\tau,i})(\tau-1)$ to $w(\hat{Q}_{\tau,i})\tau$ or vice
versa. Changing one term leads to a change of (26) of magnitude
$\frac{w(\hat{Q}_{\tau,i})}{n}$. This implies that we can get (26) within a
distance $\lvert\frac{w(\hat{Q}_{\tau,i})}{n}\lvert$ of 0. If
$w(\hat{Q}_{\tau,i})=O_{P}(1)$, which holds if $A-\hat{E}(A|L)=O_{p}(1)$ and
$\frac{1}{\hat{f}(Q_{\tau}(Y|A,L)|A,L)}>\sigma$ with probability 1 for some
$\sigma>0$, we have that (26) is $O_{P}(1/n)$. Potentially, multiple terms
change at the same time. If $M$ terms change at the same time, we have that
(26) lies within a distance $M\lvert\frac{w(\hat{Q}_{\tau,i})}{n}\lvert$ of 0.
Consequently, (26) is $O_{P}(M/n)$. If $M=o_{p}(\sqrt{n})$, we have that (26)
is $o_{P}(n^{-1/2})$.

##  Appendix C Appendix C

###  C.1 More details about the construction of estimator (19)

Assume that a (stepwise) variable selection estimator for $Q_{\tau}(Y|A,L)$
has the following form:

$\hat{Q}_{\tau}(Y|A,L)=\hat{\beta}_{\tau}A+\hat{\alpha}^{T}L^{*},$

with $L^{*}$ consisting of a subset of variables of $L$. Retargeting

$\widetilde{Q}_{\tau}(Y\mid A,L)=\hat{Q}_{\tau}(Y\mid
A,L)+\hat{\epsilon}w(\hat{Q}_{\tau,i}),$

then leads to:

$\widetilde{Q}_{\tau}(Y|A,L)=\hat{\beta}A+\hat{\alpha}^{T}L^{*}+\epsilon\frac{A-\hat{E}(A|L)}{\hat{f}(\hat{Q}_{\tau}(Y|A,L)|A,L)}.$

Then, one can obtain

$\displaystyle\widetilde{Q}_{\tau}(Y|A,L)-E\Big{(}\widetilde{Q}_{\tau}(Y|A,L)\big{\lvert}L\Big{)}$
$\displaystyle=\hat{\beta}A+\hat{\alpha}^{T}L^{*}+\epsilon\frac{A-\hat{E}(A|L)}{\hat{f}(\hat{Q}_{\tau}(Y|A,L)|A,L)}-E\left(\hat{\beta}A+\hat{\alpha}^{T}L^{*}+\epsilon\frac{A-\hat{E}(A|L)}{\hat{f}(\hat{Q}_{\tau}(Y|A,L)|A,L)}\big{\lvert}L\right)$
$\displaystyle=\hat{\beta}(A-\hat{E}(A|L))+\epsilon\left\\{\frac{A-\hat{E}(A|L)}{\hat{f}(\widetilde{Q}_{\tau}(Y|A,L)|A,L)}-E\left(\frac{A-\hat{E}(A|L)}{\hat{f}(\widetilde{Q}_{\tau}(Y|A,L)|A,L)}\right)\right\\}$

which allows us to construct a TMLE estimator as follows:

$\displaystyle\hat{\Psi}_{\tau}^{TMLE}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\frac{A_{i}-\hat{E}(A_{i}|L_{i})}{\frac{1}{n}\sum_{i=1}^{n}(A_{i}-\hat{E}(A_{i}|L_{i}))^{2}}\Bigg{[}\widetilde{Q}_{\tau}(Y_{i}|A_{i},L_{i})$
$\displaystyle\qquad\qquad-\hat{E}\Big{(}\widetilde{Q}_{\tau}(Y_{i}|A_{i},L_{i})\big{\lvert}L_{i}\Big{)}+\frac{\tau-I\\{Y_{i}\leq\widetilde{Q}_{\tau}(Y_{i}|A_{i},L_{i})\\}}{\hat{f}(\widetilde{Q}_{\tau}(Y_{i}|A_{i},L_{i})|A_{i},L_{i})}\Bigg{]}$
$\displaystyle=\hat{\beta}+\frac{1}{n}\sum_{i=1}^{n}\frac{A_{i}-\hat{E}(A_{i}|L_{i})}{\frac{1}{n}\sum_{i=1}^{n}(A_{i}-\hat{E}(A_{i}|L_{i}))^{2}}\Bigg{[}\epsilon\Bigg{\\{}\frac{A_{i}-\hat{E}(A_{i}|L_{i})}{\hat{f}(\widetilde{Q}_{\tau}(Y_{i}|A_{i},L_{i})|A_{i},L_{i})}$
$\displaystyle\qquad-\hat{E}\left(\frac{A_{i}-\hat{E}(A_{i}|L_{i})}{\hat{f}(\widetilde{Q}_{\tau}(Y_{i}|A_{i},L_{i})|A_{i},L_{i})}\big{\lvert}L_{i}\right)\Bigg{\\}}+\frac{\tau-I\\{Y_{i}\leq\widetilde{Q}_{\tau}(Y_{i}|A_{i},L_{i})\\}}{\hat{f}(\widetilde{Q}_{\tau}(Y_{i}|A_{i},L_{i})|A_{i},L_{i})}\Bigg{]}.$

Similarly to (13), the last term will be approximately zero.

##  Appendix D Appendix D

###  D.1 Additional simulation results for experiment 1

Experiment 1 was partially carried out with the computationally more expensive
‘haldensify’ estimator for the conditional density of $Y$ given $A$ and $L$.
In setting 3 we did not consider a TMLE estimator because this was
computationally too heavy. The results are shown in Table 5.

In general, the DML estimator performs better when using this conditional
density estimator. However, the TMLE estimator performs better when using the
simpler ‘FKSUM’ density estimator, especially for more extreme quantiles.

|  | $\tau=0.5$ |  | $\tau=0.75$ |  | $\tau=0.9$
---|---|---|---|---|---|---
Setting | estimator | bias | SD | SE | Cov |  | bias | SD | SE | Cov |  | bias | SD | SE | Cov
| Oracle | 0.21 | 19 | 20 | 97.4 |  | -1.2 | 33 | 35 | 95.3 |  | -2.9 | 57 | 60 | 95.7
| Plug-in | -70 | 13 | 1.5 | 0.1 |  | -73 | 14 | 1.8 | 0.1 |  | -63 | 23 | 3.5 | 0.2
| DML | -3.6 | 25 | 22 | 91.0 |  | -6.6 | 42 | 32 | 83.6 |  | -32 | 40 | 35 | 79.5
| DML-CF | -5.9 | 22 | 21 | 91.9 |  | -17 | 33 | 30 | 86.8 |  | -40 | 36 | 25 | 59.2
| TMLE | -39 | 25 | 23 | 58.8 |  | -69 | 41 | 33 | 44.7 |  | -140 | 66 | 45 | 15.1
1 | TMLE-CF | -1.4 | 24 | 24 | 95.3 |  | -1.4 | 39 | 30 | 85.6 |  | 15 | 76 | 24 | 47.8
| Oracle | 0.044 | 25 | 26 | 94.6 |  | -0.025 | 44 | 45 | 94.6 |  | -4.0 | 74 | 79 | 95.3
| Plug-in | -87 | 27 | 2.8 | 0.5 |  | -120 | 40 | 4.1 | 1.0 |  | -150 | 59 | 6.8 | 1.3
| DML | 3 | 32 | 27 | 90.7 |  | -11 | 55 | 39 | 78.5 |  | -97 | 64 | 42 | 39.3
| DML-CF | -5.0 | 28 | 24 | 88.0 |  | -34 | 43 | 33 | 74.0 |  | -110 | 60 | 26 | 15.2
| TMLE | -39 | 30 | 29 | 70.0 |  | -92 | 57 | 43 | 44.1 |  | -210 | 72 | 51 | 6.3
2 | TMLE-CF | 7.0 | 31 | 28 | 91.6 |  | 6.5 | 49 | 33 | 82.3 |  | 32 | 89 | 25 | 40.9
| Oracle | 0.19 | 3.4 | 3.6 | 94.6 |  | -0.11 | 6.2 | 6.3 | 94.0 |  | -0.30 | 10 | 11 | 94.9
| Plug-in | -17 | 6.4 | 1.6 | 1.5 |  | -25 | 8.2 | 1.9 | 0.5 |  | -39 | 11 | 2.1 | 0.3
| DML | -5.2 | 7.1 | 4.2 | 63.4 |  | -8.1 | 9.3 | 5.8 | 60.8 |  | -33 | 11 | 4.5 | 2.6
3 | DML-CF | -4.8 | 4.0 | 3.4 | 66.7 |  | -13 | 5.7 | 3.4 | 13.4 |  | -2.8 | 8.1 | 3.1 | 0.5
Table 5: Simulation results for experiment 1 with density estimated via the
‘haldensify’ R package: sample size 500, Monte Carlo bias, Monte Carlo
standard deviation (SD), average of the influence function based standard
errors (SE) and coverage of 95% Wald confidence intervals (Cov). All values
have been multiplied by $10^{2}$.

Table 6 shows a summary of the term we tried to make zero with targeted
learning (11). When using cross-fitting, we observe that this term is very
close to zero, averaging over all simulation runs.

|  | TMLE | TMLE-CF
---|---|---|---
Minimum |  | -0.22 | -0.059
1st Quartile |  | -0.087 | -0.0029
Median |  | -0.0028 | 0.00013
Mean |  | -0.036 | 0.00023
3rd Quartile |  | 0.00013 | 0.0033
Maximum |  | 0.20 | 0.064
Table 6: Summary of the term we wish to make zero with targeted learning (11) for the homoscedastic setting in experiment 1. Quantile | estimator |  | bias | SD | SE | Cov
---|---|---|---|---|---|---
$0.5$ | TMLE.1 |  | -0.38 | 0.25 | 0.19 | 44.6
TMLE.1-CF |  | 0.082 | 0.27 | 0.18 | 83.0
$0.75$ | TMLE.1 |  | -0.71 | 0.34 | 0.25 | 25.2
TMLE.1-CF |  | 0.10 | 0.44 | 0.25 | 73.7
$0.9$ | TMLE.1 |  | -1.1 | 0.47 | 0.27 | 14.0
TMLE.1-CF |  | 0.16 | 0.75 | 0.34 | 62.2
Table 7: Simulation results for the homoscedasticity setting in experiment 1
with the 1-step TMLE estimator: sample size $n=500$, quantile $\tau$, Monte
Carlo bias, Monte Carlo standard deviation (SD), average of the influence
function based standard errors (SE) and coverage of 95% Wald confidence
intervals (Cov).

The TMLE(-CF) estimator for continuous exposure was adjusted such that only 1
targeting step was carried out. For binary exposure, this step was repeated
until convergence. In Table 7, we show the results of the 1-step TMLE
estimator for binary exposure (experiment 1, setting 1). Comparing these
results with the results in Table 1, it can be observed that 1-step targeting
does not achieve the same performance as targeting until convergence.

###  D.2 Additional simulation results for experiment 4

Table 8 contains the results for the following estimators for the data
generating mechanism described in experiment 4: the oracle estimator (Oracle),
the naive plug-in (Plug-in) estimator, the debiased machine learning estimator
(8) based on quantile random forests without cross-fitting (DML) and with
5-fold cross-fitting (DML-CF), the targeted learning estimator (13) without
cross-fitting (TMLE) and with 5-fold cross-fitting (TMLE-CF). We furthermore
included some adjusted estimators. One is obtained by performing quantile
regression with AIC-based backward stepwise variable selection (QRvs). Another
one is the debiased machine learning estimator defined in (18) where quantile
random forests are replaced by quantile regression with AIC-based backward
stepwise variable selection, once without cross-fitting (DMLvs) and once with
5-fold cross-fitting (DMLvs-CF). The third one is the TMLE estimator defined
in (19), once without cross-fitting (TMLEvs) and once with 5-fold cross-
fitting (TMLEvs-CF).

|  |  | $\mathbf{n=250}$ |  | $\mathbf{n=500}$
---|---|---|---|---|---
Quantile | estimator |  | bias | SD | SE | Cov |  | bias | SD | SE | Cov
$\tau=0.5$ | Oracle |  | 0.66 | 15 | 16 | 96.0 |  | 0.25 | 10 | 11 | 95.7
QRvs |  | 2.4 | 17 | 16 | 91.9 |  | 0.36 | 12 | 11 | 92.8
Plug-in |  | -36 | 11 | 3.6 | 0.9 |  | -23 | 8.6 | 2.6 | 2.7
DML |  | -12 | 20 | 20 | 90.1 |  | -20 | 13 | 12 | 59.4
DML-CF |  | 6.6 | 16 | 15 | 90.9 |  | 4.0 | 10 | 10 | 92.1
DMLvs |  | 2.0 | 18 | 13 | 81.8 |  | 0.33 | 12 | 9.9 | 90.2
DMLvs-CF |  | 1.6 | 17 | 9.8 | 74.0 |  | 0.043 | 12 | 8.5 | 84.1
TMLE |  | -33 | 18 | 20 | 61.0 |  | -26 | 15 | 12 | 40.0
TMLE-CF |  | 19 | 20 | 22 | 87.7 |  | 8.5 | 13 | 15 | 92.8
TMLEvs |  | -29 | 25 | 13 | 67.5 |  | -0.49 | 15 | 10 | 81.7
TMLEvs-CF |  | 1.3 | 20 | 19 | 94.3 |  | -0.32 | 13 | 12 | 94.8
$\tau=0.75$ | Oracle |  | -0.016 | 16 | 18 | 96.7 |  | -0.54 | 12 | 12 | 95.0
QRvs |  | -1.4 | 19 | 17 | 91.3 |  | 0.0095 | 13 | 12 | 93.2
Plug-in |  | -42 | 12 | 4.1 | 1.1 |  | -29 | 9.1 | 2.9 | 1.4
DML |  | -19 | 20 | 19 | 77.8 |  | -22 | 13 | 11 | 47.3
DML-CF |  | -0.44 | 18 | 17 | 94.4 |  | -2.3 | 12 | 11 | 92.3
DMLvs |  | -0.78 | 20 | 12 | 77.2 |  | -0.24 | 13 | 9.9 | 86.2
DMLvs-CF |  | -0.90 | 18 | 10 | 72.7 |  | -0.21 | 12 | 8.8 | 82.5
TMLE |  | -41 | 21 | 20 | 45.0 |  | -37 | 16 | 12 | 19.9
TMLE-CF |  | 16 | 24 | 24 | 88.5 |  | 5.7 | 16 | 16 | 93.7
TMLEvs |  | -3.9 | 28 | 13 | 57.4 |  | -1.1 | 17 | 9.9 | 72.5
TMLEvs-CF |  | 0.52 | 21 | 21 | 94.0 |  | -0.28 | 14 | 13 | 93.5
$\tau=0.9$ | Oracle |  | -0.88 | 21 | 22 | 95.7 |  | 0.39 | 14 | 15 | 95.8
QRvs |  | 0.065 | 23 | 20 | 89.7 |  | 0.094 | 16 | 14 | 90.8
Plug-in |  | -64 | 8.4 | 037 | 0.0 |  | -50 | 7.5 | 3.1 | 0.0
DML |  | -50 | 16 | 19 | 27.3 |  | -46 | 9.1 | 9.1 | 0.8
DML-CF |  | -9.0 | 25 | 26 | 93.0 |  | -14 | 16 | 16 | 84.8
DMLvs |  | -0.52 | 25 | 11 | 60.9 |  | -0.23 | 17 | 9.3 | 73.8
DMLvs-CF |  | 0.41 | 20 | 11 | 73.0 |  | 0.021 | 14 | 9.5 | 80.7
TMLE |  | -67 | 23 | 23 | 18.1 |  | -67 | 16 | 13 | 3.1
TMLE-CF |  | 21 | 38 | 33 | 83.1 |  | 4.2 | 27 | 22 | 88.1
TMLEvs |  | -6.1 | 41 | 11 | 30.4 |  | -1.2 | 26 | 9.5 | 42.8
TMLEvs-CF |  | 0.21 | 24 | 24 | 94.4 |  | -0.18 | 16 | 16 | 93.7
Table 8: Simulation results for variable selection: sample size $n$, quantile
$\tau$, Monte Carlo bias, Monte Carlo standard deviation (SD), average of the
influence function based standard errors (SE) and coverage of 95% Wald
confidence intervals (Cov). All values have been multiplied by $10^{2}$.

###  D.3 Propensity score comparison: experiment 1 vs. 2

To get some insight into the propensity scores used in experiment 2,
histograms are provided below. The left histogram shows a histogram of the
propensity scores used in experiment 1 (with binary $A$). On the right, a
histogram of the propensity scores from experiment 2 is shown. Each time, a
sample of size 100,000 was generated according to the two data-generating
mechanisms. Then a histogram of the propensity scores was constructed. It is
clear that more extreme propensity scores occur in experiment 2.

Figure 2: Propensity score histograms corresponding to the binary exposure
settings of experiment 1 (left) and experiment 2 (right).

##  Appendix E Appendix E

###  E.1 Sensitivity analysis for the data analysis

As shown by the simulation studies, the uncertainty of the proposed TMLE
estimator is sometimes underestimated for extreme quantiles. Depending on how
the folds in the cross-fitting procedure are chosen, the results may differ.
We repeated the analysis for the 90th percentile 10 times with different
seeds, once with 5-fold cross-fitting and once with 10-fold cross-fitting. The
results are shown in Fig. 3

Figure 3: Repeating the analysis 10 times for the 90th percentile, once with 5
folds (shown in black), once with 10 folds (shown in brown). The blue interval
corresponds to the one reported in Section 5.

##  Appendix F Appendix F

###  F.1 Generalization to link function

There are settings in which it could be useful to generalize the estimand (3)
(or, equivalently, (2) so that it includes a link function. A realistic
example is a survival analysis where the outcome can only take positive values
(e.g. when it is assumed to be exponentially distributed with exponential
hazard). The quantile in our estimand could take negative values, which makes
no sense in this setting. It would be reasonable here to use a link function
and replace model (1) by

$g\\{Q_{\tau}(Y|A,L)\\}=\beta_{\tau}A+\omega_{\tau}(L),$ (27)

with link function $g$ equal to the natural logarithm. This straightforwardly
leads to the following estimand for $\beta_{\tau}$:

$\Theta_{\tau}=\frac{E[\\{A-E(A|L)\\}(g\\{Q_{\tau}(Y|A,L)\\}-E[g\\{Q_{\tau}(Y|A,L)\\}|L])]}{E[\\{A-E(A|L)\\}^{2}]}.$
(28)

For this estimand, the efficient influence function can be calculated in the
same way as in the proof of Theorem 1. The main difference occurs in Lemma 2.
With link function, by making use of the chain rule one would obtain

$\displaystyle\frac{d}{dt}(A-E_{t}(A|L))g\\{Q_{\tau,t}(Y|A,L)\\}\Big{\lvert}_{t=0}$
$\displaystyle=\frac{d}{dt}(A-E_{t}(A|L))\Big{\lvert}_{t=0}g\\{Q_{\tau,t}(Y|A,L)\\}+(A-E(A|L))\frac{d}{dt}\left(g\\{Q_{\tau,t}(Y|A,L)\\}\right)\Big{\lvert}_{t=0}$
$\displaystyle=-\frac{\mathbbm{1}_{\tilde{l}}(L)}{f(L)}(\tilde{a}-E(A|L))g\\{Q_{\tau,t}(Y|A,L)\\}+\\{A-E(A|L)\\}g^{\prime}\\{Q_{\tau}(Y|A,L)\\}\frac{d}{dt}Q_{\tau,t}(Y|A,L)\Big{\lvert}_{t=0}$
$\displaystyle=-\frac{\mathbbm{1}_{\tilde{l}}(L)}{f(L)}(\tilde{a}-E(A|L))g\\{Q_{\tau,t}(Y|A,L)\\}+\\{A-E(A|L)\\}g^{\prime}\\{Q_{\tau}(Y|A,L)\\}\mathbbm{1}_{\tilde{a},\tilde{l}}(A,L)\frac{\tau-I\\{\tilde{y}\leq
Q_{\tau}(Y|A,L)\\}}{f(Q_{\tau}(Y|A,L),A,L)}$

such that the efficient influence function is

$\displaystyle\phi(Y,A,L)=\frac{A-E(A|L)}{E[\\{A-E(A|L)\\}^{2}]}\bigg{[}g\\{Q_{\tau}(Y|A,L)\\}-E[(g\\{Q_{\tau}(Y|A,L)\\}|L]$
$\displaystyle\qquad\qquad\qquad+\frac{\tau-I\\{Y\leq
Q_{\tau}(Y|A,L)\\}}{f\\{Q_{\tau}(Y|A,L)|A,L\\}}g^{\prime}\\{Q_{\tau}(Y|A,L)\\}-\Theta_{\tau}\\{A-E(A|L)\\}\bigg{]}.$

The proposed estimators can be easily adjusted for the estimation of
$\Theta_{\tau}$.