# Independence weights for causal inference with continuous exposures

Jared D. Huling1, Noah Greifer2, Guanhua Chen3

1Division of Biostatistics, University of Minnesota
2Department of Mental Health,
Johns Hopkins Bloomberg School of Public Health
3Department of Biostatistics and Medical Informatics,
University of Wisconsin-Madison
corresponding authors<EMAIL_ADDRESS>authors<EMAIL_ADDRESS>
###### Abstract

Studying causal effects of continuous exposures is important for gaining a
deeper understanding of many interventions, policies, or medications, yet
researchers are often left with observational studies for doing so. In the
observational setting, confounding is a barrier to estimation of causal
effects. Weighting approaches seek to control for confounding by reweighting
samples so that confounders are comparable across different values of the
exposure, yet for continuous exposures, weighting methods are highly sensitive
to model misspecification. In this paper we elucidate the key property that
makes weights effective in estimating causal quantities involving continuous
exposures. We show that to eliminate confounding, weights should make exposure
and confounders independent on the weighted scale. We develop a measure that
characterizes the degree to which a set of weights induces such independence.
Further, we propose a new model-free method for weight estimation by
optimizing our measure. We study the theoretical properties of our measure and
our weights, and prove that our weights can explicitly mitigate exposure-
confounder dependence. The empirical effectiveness of our approach is
demonstrated in a suite of challenging numerical experiments, where we find
that our weights are quite robust and work well under a broad range of
settings.
Keywords— Continuous treatments, Distance Covariance, Confounding, Covariate
Balance, Balancing Weights

##  1 Introduction

Confounding is a major barrier to studying causal effects of exposures from
observational data, regardless of whether the exposure is binary, continuous,
or otherwise. Considerable work has focused on the development of approaches
for studying causal effects of binary or otherwise discrete-valued exposures
from observational data. With continuous exposures, however, the choices of
methods for confounding control are far more limited, and clear guidance that
can help practitioners choose among the available methods is lacking. A common
approach to reduce confounding by observed variables is using the propensity
score, which was initially proposed for binary treatments Rosenbaum and Rubin
(1983) and has been generalized to the setting of continuous treatments
(Hirano and Imbens, 2004; Imai and Van Dyk, 2004; Zhu et al., 2015; Galvao and
Wang, 2015). With binary treatments, the causal effect of interest is often
the average treatment effect, which can be estimated as a difference in the
weighted averages of the treatment group outcomes, where the weights are
proportional to the inverse of the propensity score; this method is known as
inverse probability weighting (IPW). With continuous treatments, the interest
is often in estimation of the causal dose-response functionals (Robins et al.,
2000; van der Laan and Robins, 2003) such as the the causal average dose-
response function (ADRF), which can be estimated using a weighted regression
of the outcome on the treatment, where the weights are proportional to the
inverse of the conditional density of the treatment given the covariates, the
generalized propensity score (GPS).

In the binary treatment setting, IPW estimators can be unstable due to extreme
weights and susceptible to model misspecification (Kang and Schafer, 2007).
These issues carry over to IPW estimators for continuous treatments and are
substantially more challenging to address. The reason is that IPW estimation
via the GPS requires inverse weighting by a conditional density estimate, not
just a conditional probability. Even if the conditional mean of the treatment
given covariates is correctly specified, GPS weights can fail to perform well
if the distribution of the conditional density is misspecified (Naimi et al.,
2014). The difficulty of correctly specifying or accurately estimating a full
conditional distribution is exacerbated with increased dimension of the
confounders. The tailored use of flexible machine learning estimation
approaches such as that proposed in Zhu et al. (2015) can in many cases yield
substantial improvements, however, as shown in our simulation studies and data
analysis, they are not fully immune from the issues of GPS weighting and can
perform poorly and/or yield large weights in practice.

As misspecification of a conditional density model can be difficult to
diagnose and assess, several works have focused on directly estimating weights
through various optimization criteria with the goal of inducing a lack of
correlation between marginal moments of confounders and the exposure.
Approaches along this line of work include the generalized covariate balancing
propensity score (CBPS) approach of Fong et al. (2018), covariate association
eliminating weights (Yiu and Su, 2018), and entropy balancing weights
(Tübbicke, 2020; Vegetabile et al., 2021). Because these approaches focus on
estimation of weights directly as opposed to estimating a conditional density
explicitly and inverting it, they tend to be more robust and effective
empirically than direct modeling of the GPS. While intuitively appealing,
these approaches require careful choices of which moments of both the
confounders and the exposure to “decorrelate”. Yet, there is no guidance on
how to specify the correct moments necessary to mitigate bias in estimation of
the ADRF due to confounding. Missing important moments can leave substantial
residual dependence between confounder and exposure, see e.g. Figure 1. In our
simulations, we show that the choice of moments is indeed critically important
in practice and that numerical instability can arise when too many moments are
used. The tension between including enough moments to reduce bias and the
instability of weights as more moments are included can make the use of these
methods challenging in practice. Kernel Optimal Orthogonality Weighting (KOW)
(Kallus and Santacatterina, 2019) and a generalization of KOW proposed in
Martinet (2020), are kernel-based nonparametric extensions of direct weights
estimation ideas. They focus on estimating weights to decorrelate over
function spaces of exposure and confounders such that there is no need of
choosing models or what moments to decorrelate. Yet, careful tuning is still
required when flexible kernels are used, and, on the other hand, when
inflexible kernels are used, there is no guarantee that the resulting weights
fully mitigate confounding. Thus, the user is often left with a difficult
choice of which kernel to use and unclear guidance on how to choose the
kernel’s tuning parameters.

Our work aims to achieve several goals. First, we provide clarity on the role
of weights in estimation of the ADRF. To do so, we provide a general
decomposition of the error of a weighted nonparametric estimator of the ADRF
and demonstrate that, under broad conditions, the ideal weights should induce
independence between the treatment variable and confounders to guarantee
mitigation of confounding bias when estimating the ADRF. While this is already
intuitively understood in the literature, our decomposition shows the precise
form that dependence shows up in the error of an estimator. We also show that
dependence plays a key role in other estimands, such as the causal quantile
dose-response function. Second, we develop a measure based on energy
statistics (Székely et al., 2007; Székely and Rizzo, 2013) that allows one to
assess how well a set of weights is able to induce independence, where smaller
values of our measure indicate the weights yield less treatment-confounder
dependence and a value of zero indicates complete independence between the
treatment and confounders in the weighted data. Our proposed measure can
effectively characterize how well dependence is mitigated by a set of weights
without the need for fine-tuning hyperparameters. Finally, we propose a new
approach for estimating weights, which we call the distance covariance optimal
weights (DCOWs), by optimizing our measure. Thus, the proposed weights
directly aim to mitigate dependence, the key source of error in a weighted
nonparametric estimate of the ADRF. Our weight construction approach does not
require modeling a conditional density, careful tuning of hyperparameters, or
choosing which moments of covariates and treatment to decorrelate, which means
it can be readily and effectively used “out-of-the-box” by practitioners with
varying degrees of statistical sophistication. Because our approach targets
the key source of error in estimates of causal estimands involving continuous
exposures, our weights yield estimators of the ADRF that are highly robust
empirically and work well in a variety of settings, including when confounders
are of moderate dimension. The key benefits of our approach are its ease of
use, requiring minimal to no fine-tuning on the part of the user, and its
ability to provide highly reliable estimates in many settings. As such, it has
the potential for wide applicability in practice.

We provide some theoretical results for our proposal, showing that our weights
indeed induce independence between treatment and covariates and that, with a
small penalty on the variability of the weights, our proposal results in the
same convergence rate as a nonparametric regression estimate of the ADRF in a
scenario with no confounding. Although adding a penalty to reduce weight
variability involves the inclusion of a tuning parameter, our proposed
approach rarely results in weights with large variability even without
penalization, so careful tuning of the parameter that controls weight
variability is rarely necessary, as evidenced by our simulation studies, which
investigate a wide variety of scenarios with strong and complex confounding
and scenarios with moderately high-dimensional confounding, in all of which we
allow the tuning parameter to remain at its default value.

Our proposed weights can be used beyond simple weighted nonparametric
estimators of the ADRF. Recently, Kennedy et al. (2017) and Díaz and van der
Laan (2013) extended the idea of doubly robust estimation to continuous
treatments, allowing for estimates that combine outcome regression models
(Imbens, 2004; Hill, 2011) and the conditional density models. This allows for
relaxed dependence on the correctness of the regression and conditional
density models. However, doubly robust estimators are not immune to highly
variable weights and their finite sample performance can suffer materially if
the conditional density model is misspecified. We show that our weights can
enhance doubly-robust estimators. In particular, our measure of dependence can
provide a sharper bound on error for doubly robust estimators. This reveals
that pairing of a reasonable outcome regression model with our weights in a
doubly-robust fashion can be effective for estimating the ADRF. Our weight
estimation procedure can also be modified to allow for exact decorrelation of
chosen moments, allowing one to mitigate dependence while emphasizing
decorrelation of certain features known to be important. This may improve
performance in some settings.

The remainder of this paper is organized as follows. We investigate the role
of dependence in weighted estimation of the ADRF in Section 2 and we develop a
criterion that assesses how much dependence is mitigated by a set of weights,
propose a new weight estimation strategy, and provide some corresponding
theory in Section 3. We demonstrate the effectiveness of our approach in
finite samples with a suite of challenging simulation studies in Section 4 and
illustrate the use of our approach in a real-world study of EHR data in
Section 5. We conclude the paper with some discussion.

##  2 Continuous Treatments and Confounding, Weighting and Dependence

###  2.1 Setup, notation, and assumptions

The observable quantities we consider consist of the random triplet $({\bf
X},A,Y)$, where ${\bf X}\in\mathcal{X}\subseteq\mathbb{R}^{p}$ is a vector of
pre-treatment covariates, $A\in\mathcal{A}\subseteq\mathbb{R}$ is a
continuous-valued treatment variable indicating the assigned dose for a unit,
and $Y\in\mathcal{Y}\subseteq\mathbb{R}$ is an outcome of interest. The
variate $({\bf X},A,Y)$ has a joint distribution $F_{{\bf X},A,Y}$ with
respect to a dominating measure. We denote the marginal density of the
treatment and covariates as $f_{A}(a)$ and $f_{\bf X}({\bf x})$, respectively,
the conditional density of the treatment given ${\bf X}$ as $f_{A|{\bf
X}}(a|{\bf x})$, and their joint density as $f_{{\bf X},A}(a,{\bf x})$.
Similarly, corresponding distribution functions are denoted
$F_{A}(a)=\mathbb{P}(A\leq a)$, $F_{\bf X}({\bf x})=\mathbb{P}({\bf X}\leq{\bf
x})$, $F_{A|{\bf X}}(a|{\bf x})=\mathbb{P}(A\leq a\>|\>{\bf X}={\bf x})$, and
$F_{{\bf X},A}({\bf x},a)=\mathbb{P}({\bf X}\leq{\bf x},A\leq a)$. Our
observed data consists of $n$ i.i.d. samples $({\bf
X}_{i},A_{i},Y_{i})_{i=1}^{n}$ from $({\bf X},A,Y)$. Note that we drop the
subscripts on the density and cumulative distribution functions when there is
no ambiguity.

To clarify the causal quantity of interest, we work under the potential
outcomes framework, wherein the potential outcome function $Y(a)$ for
$a\in\mathcal{A}$ is the outcome that would be observed if $A$ were set to the
value $a$. The causal quantity of prime interest in this paper is the mean
potential outcome function, also called the causal average dose-response
function (ADRF), which is

$\mu(a)\equiv\mathbb{E}\left[Y(a)\right],\text{ for }a\in\mathcal{A}.$

The causal ADRF $\mu(a)$ can be identified, or expressed in terms of
observational data, under standard causal assumptions. These assumptions,
which we employ throughout this paper are 1) consistency, which posits that
$A=a$ implies $Y=Y(a)$, 2) positivity, which states that all values of the
treatment as possible across the covariate space in that $f(a|{\bf X}={\bf
x})\geq\nu>0$ for all ${\bf x}\in\mathcal{X}$ for some constant $\nu$, and 3)
ignorability of the assignment mechanism: $Y(a){\perp\\!\\!\\!\perp}A\>|\>{\bf
X}$ for all $a\in\mathcal{A}$, where ${\perp\\!\\!\\!\perp}$ denotes
(conditional) independence (Dawid, 1979). Under assumptions 1)-3), the dose-
response function is identified as

$\mu(a)=\mathbb{E}_{\bf X}\left(\mathbb{E}\left[Y\>|\>{\bf
X},A=a\right]\right)=\mathbb{E}_{\bf X}\left[\mu({\bf X},a)\right].$

Estimation of the dose-response function via regression-based estimation of
the mean function $\mu({\bf X},a)\equiv\mathbb{E}\left[Y\>|\>{\bf
X},A=a\right]$ can be highly challenging. Misspecification of the regression
function can result in poor estimation of $\mu(a)$, and nonparametric
estimation of the regression function is also difficult, especially when ${\bf
X}$ is not low dimensional. Instead, in this paper we focus on weighting-based
estimators of $\mu(a)$. A benefit of weighting estimators is that the dose-
response function can be flexibly estimated by univariate (weighted)
nonparametric regression, whereas with regression-based estimation, the need
to incorporate confounders in a regression may make flexible estimation of the
ADRF more difficult due to the additional dimensions. A conceptual benefit of
weighting methods over regression-based approaches is that they provide a
clear separation between the design phase and the analysis phase of a study.
This separation is critical when substantial or iterative modeling effort is
required to control for confounding.

The focus of this manuscript is in developing a criterion that can evaluate
how effective a set of weights is in mitigating confounding and providing a
new set of weights which is optimal under the criterion. We will demonstrate
that this approach results in weights that yield robust control of the
dependence between $A$ and ${\bf X}$ and thus robust estimation of the ADRF.

###  2.2 Ignorability, independence, and the GPS

To reliably estimate the causal dose response function using observational
data, sources of structural bias should be mitigated, among which the
confounding bias is the most common. Consider the following DAG which
illustrates our setting, where the covariate ${\bf X}$ is a confounder in
studying the causal relationship between $A$ and $Y$. To mitigate the
confounding, the backdoor path $A\leftarrow{\bf X}\rightarrow Y$ must be
blocked. One way of blocking the backdoor path is by removing the arrow
between $A$ and ${\bf X}$ (i.e. making $A$ independent of ${\bf X}$). In a
randomized trial, the independence between $A$ and ${\bf X}$ holds due to
randomization. This motivates us to create a pseudo-population mimicking the
one we would observe under such a trial by reweighting the subjects in an
observational study such that $A$ is independent of ${\bf X}$ in the pseudo-
population. This is equivalent to finding a function of ${\bf X}$, $r({\bf
X})$, such that $A{\perp\\!\\!\\!\perp}{\bf X}\>|\>r({\bf X})$. Together with
the ignorability condition, if $A{\perp\\!\\!\\!\perp}{\bf X}\>|\>r({\bf X})$,
$r({\bf X})$ is also a balancing score such that weak unconfoundness holds
given the balancing score, i.e., for every $a$, $f_{A}(a\>|\>r(a,{\bf
X}),Y(a))=f_{A}(a\>|\>r(a,{\bf X}))$ (Hirano and Imbens, 2004). The
generalized propensity score $f_{A|{\bf X}}(A\>|\>{\bf X})$ (GPS) (Hirano and
Imbens, 2004) is such a balancing score and extends the pioneering work of
Rosenbaum and Rubin (1983) proposed for binary treatments to continuous
treatments.

For categorical treatments, weights for inducing independence between $A$ and
${\bf X}$ are computed as $1/\mathbb{P}(A=a\>|\>{\bf X})$, with the
variability-reducing stabilized weights computed as
$\mathbb{P}(A=a)/\mathbb{P}(A=a\>|\>{\bf X})$. For continuous treatments,
stabilized GPS weights are computed as $f_{A}(A)/f_{A|{\bf X}}(A\>|\>{\bf X})$
(Robins et al., 2000), which naturally arise in estimating equations for the
dose-response function via semiparametric theory (Kennedy et al., 2017).
Estimating the weights requires correct specification not only of the mean of
the conditional density of the treatment, but also of its shape. When either
of these are misspecified, bias can result in estimates of the ADRF (Naimi et
al., 2014; Zhu et al., 2015), indicating a particularly sensitive reliance on
correct modeling of the full relationship between covariates and the
conditional distribution of the treatment. Furthermore, similar to standard
propensity score weights, GPS modeling can yield extreme weights, leading to
unstable estimation. Weight trimming/capping may alleviate the problem of
large weights but can be seen as ad hoc and may change the estimand (Crump et
al., 2009): the estimated ADRF will correspond to the population represented
by the newly weighted sample rather than to the original target population.

The stabilized GPS weights $f_{A}(A)/f_{A|{\bf X}}(A\>|\>{\bf X})$ have
several key properties that have motivated work to improve upon the GPS
weights. Namely, when weighting by the GPS weights in the population sense,
they 1) result in independence of $A$ and ${\bf X}$, 2) preserve the marginal
distributions of ${\bf X}$ and $A$, and 3) have mean 1. These properties are
listed more explicitly in Appendix A. Instead of indirectly estimating the
weights by estimating the GPS, a nascent line of work has involved methods
which directly estimate weights designed to satisfy the above three properties
via balancing criteria. For example, work has focused on estimation of weights
that induce zero marginal correlation between treatment and covariates subject
to constraints that weighted marginal means of covariates and treatment are
the same as the unweighted marginal means (Fong et al., 2018; Yiu and Su,
2018; Vegetabile et al., 2021). Although they are more robust than GPS
weights, these approaches rely on both the correct choice of moments of the
confounders and the choice of moments of the exposure variable to decorrelate.
Yet, there is little guidance in deciding which set of moments in both
confounders and exposure to focus on.

The aforementioned three properties of GPS weights are intuitively appealing,
but it is not immediately clear which of or in what manner these properties
are important in mitigating bias in a weighted nonparametric estimator of the
ADRF. To justify which properties are crucial in weights estimation, we derive
the relationship between the properties of weighted dose-response function
estimators using generic weights and the systematic source of error of the
weighted estimator for the ADRF. We demonstrate that the ability of a set of
weights to induce independence between $A$ and ${\bf X}$ is critical for
reducing the bias of a weighted estimator. Thus, we show that inducing
independence should be the primary aim when estimating weights to study the
causal effect of a continuous exposure, i.e. treatment.

###  2.3 A general error decomposition for weighted nonparametric estimators
of the ADRF

In this section we aim to provide an explicit mechanistic connection between
dependence and bias in weighted estimates of the ADRF. Although it is
understood that using weights constructed from a well-estimated conditional
density are consistent (Kennedy et al., 2017), it is unclear what role weights
play more generally. It is not immediately clear what the connection is
between the balance criteria that aim to “decorrelate” moments of covariates
and moments of the treatment variable and the systematic bias in estimating
the ADRF. We will investigate the precise source of the systematic bias of an
estimator and illuminate the role of the weights in influencing the bias.
Here, we focus on weighted Nadaraya-Watson estimators of the ADRF for clarity
of presentation, though the key message applies to weighted local polynomial
regression and other weighted nonparametric regression.

The response can be expressed as $Y_{i}=\mu({\bf
X}_{i},A_{i})+\varepsilon_{i}$, where $\varepsilon_{i}\equiv
Y_{i}(A_{i})-\mu({\bf X}_{i},A_{i})$. By construction, $\varepsilon_{i}$ have
mean zero but are not necessarily identically distributed. Given any set of
weights ${\bf w}=(w_{1},\dots,w_{n})$ and a kernel
$K_{h}(A_{i}-a_{0})=K(\frac{A_{i}-a_{0}}{h})/h$ centered at $A=a_{0}$ with
bandwidth $h>0$, the weighted Nadaraya-Watson estimator of the dose-response
function at $A=a_{0}$ is

$\mu_{NW}^{{\bf
w}}(a_{0})=\frac{\sum_{i=1}^{n}Y_{i}w_{i}K_{h}(A_{i}-a_{0})}{\sum_{i=1}^{n}K_{h}(A_{i}-a_{0})}.$
(1)

Given any ${\bf w}$, the error of (1) at $A=a_{0}$ can be decomposed as

$\displaystyle\widehat{\mu}_{NW}^{{\bf w}}(a_{0})-\mu(a_{0})={}$
$\displaystyle\int_{\mathcal{X}}\int_{\mathcal{A}}\mu({\bf
x},a_{0})\mathrm{d}\left[F^{n}_{{\bf X},A,{\bf w}}-F^{n}_{{\bf
X}}F^{n}_{A}\right]({\bf x},a)$ (2) $\displaystyle+\int_{\mathcal{X}}\mu({\bf
x},a_{0})\mathrm{d}\left[F^{n}_{{\bf X}}-F_{\bf X}\right]({\bf
x})+\left(\frac{f_{A}(a_{0})}{\widehat{f}^{n}_{A,h}(a_{0})}-1\right)\int_{\mathcal{X}}\mu({\bf
x},a_{0})\mathrm{d}F^{n}_{{\bf X}}({\bf x})$
$\displaystyle+\left(\frac{f_{A}(a_{0})}{\widehat{f}^{n}_{A,h}(a_{0})}-1\right)\int_{\mathcal{X}}\int_{\mathcal{A}}\mu({\bf
x},a_{0})\mathrm{d}\left[F^{n}_{{\bf X},A,{\bf w}}-F^{n}_{{\bf
X}}F^{n}_{A}\right]({\bf x},a)$
$\displaystyle+{\widehat{f}^{n}{}}^{-1}_{A,h}(a_{0})\int_{\mathcal{X}}\int_{\mathcal{A}}\left[\mu({\bf
x},a)K_{h}(a-a_{0})-\mu({\bf x},a_{0})f_{A}(a_{0})\right]\mathrm{d}F^{n}_{{\bf
X},A,{\bf w}}({\bf x},a)$
$\displaystyle+\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}w_{i}{\widehat{f}^{n}{}}^{-1}_{A,h}(a_{0}){K_{h}(A_{i}-a_{0})},$
(3)

where
$\widehat{f}^{n}_{A,h}(a_{0})=\int_{\mathcal{A}}K_{h}(a-a_{0})\mathrm{d}F^{n}_{A}(a)$
is a kernel density estimate of $f_{A}(a_{0})$. We provide a derivation of
this decomposition in the Supplementary Material. The second term on the right
is due to sampling variability only and has mean zero and converges at rate
$n^{-1/2}$ if the sample is representative of the super-population. The
expectations of the third, fourth, and fifth terms on the right above go to 0
when $h\to 0$ regardless of the weights ${\bf w}$. The last term has mean zero
regardless of both the weights and the bandwidth, though its variability is
impacted by the weights. A subtle point here is that the error decomposition
above holds if $F^{n}_{A}$ in (2) is replaced with any valid marginal
distribution function of the treatment defined on $\mathbb{R}$.

On the other hand, the first term (2) on the right above is the source of
systematic bias of the weighted estimator unrelated to kernel smoothing. In
other words, taking limits of the bandwidth of the kernel to 0 and sample size
to infinity does not make (2) vanish. This term also provides insight into why
using the weights we propose later perform well in finite sample settings when
used in treatment effect estimators, as targeting this term specifically can
help decrease the magnitude of the systematic component of the finite sample
error of an estimator. If a given set of weights induces finite-sample
independence of ${\bf X}$ and $A$ in the sense that $F^{n}_{{\bf X},A,{\bf
w}}({\bf x},a)=F^{n}_{{\bf X}}({\bf x})F^{n}_{A}(a)$ for all ${\bf
x}\in\mathcal{X}$, $a\in\mathcal{A}$, then the systematic source of error of
$\hat{\mu}^{{\bf w}}(a_{0})$ will be zero. The mean-squared error of the
estimator will, however, depend primarily on both the bias term (2) and the
variance of (3). Mitigating the variance of (3) merely amounts to controlling
the squares of the weights; however, providing a measure that can characterize
(2) is non-trivial and none exists in the literature. The systematic component
(2) of the finite sample error is bounded by the distance between $F^{n}_{{\bf
X},A,{\bf w}}({\bf x},a)$ and $F^{n}_{{\bf X}}({\bf x})F^{n}_{A}(a)$ provided
that $\mu({\bf x},a_{0})$ is bounded. The mean function $\mu({\bf x},a_{0})$
is unknowable, so without modeling the response function, constructing a
measure that bounds (2) is critical for assessing a set of weights.

In practice, the estimator $\widehat{\mu}_{NW}^{{\bf w}}(a_{0})$ in (1) may be
unstable, as the weights only appear in the numerator, so the estimated ADRF
may lie outside the range of the observed values of the response. Instead, a
more stable estimator is the following weighted average of the responses

$\mu_{NWs}^{{\bf
w}}(a_{0})=\frac{\sum_{i=1}^{n}Y_{i}w_{i}K_{h}(A_{i}-a_{0})}{\sum_{i=1}^{n}w_{i}K_{h}(A_{i}-a_{0})},$
(4)

which can be viewed as the minimizer of a weighted least squares problem where
the $i$th weight is $w_{i}K_{h}(A_{i}-a_{0})$. The estimator (4) is also a
valid estimator of the ADRF as long as the denominator is a consistent
estimator of $f_{A}(a_{0})$; in this case the key source of systematic bias
still depends on the term (2).

###  2.4 Error decompositions for estimands beyond the ADRF

Although we do not present a general theory, we emphasize that the role of
weights in their ability to induce independence between treatment and
covariates is not unique to estimation of the ADRF and applies to a wide
variety of estimands. Consider estimation of the causal dose-response quantile
function $q_{Y(a_{0})}(\alpha)=\inf\\{y:F_{Y(a_{0})}(y)\leq\alpha\\}$, where
$F_{Y(a_{0})}(y)=\mathbb{P}(Y(a_{0})\leq y)=\mathbb{E}_{{\bf
X}}\left\\{\mathbb{P}(Y\leq y|{\bf X},A=a_{0})\right\\}=\mathbb{E}_{{\bf
X}}\left\\{F_{Y|{\bf X},A}(y|{\bf X},A=a_{0})\right\\}$.

The weighted nonparametric estimator of the CDF of $Y(a_{0})$ is

$\displaystyle\widehat{F}_{Y(a_{0})}^{{\bf w}}(y)={}$
$\displaystyle\frac{\sum_{i=1}^{n}I(Y_{i}\leq
y)w_{i}K_{h}(A_{i}-a_{0})}{\sum_{i=1}^{n}K_{h}(A_{i}-a_{0})}$ (5)

and the resulting kernel quantile estimator of
$q_{Y(a_{0})}(\alpha)=\inf\\{y:\alpha\leq F_{Y(a_{0})}(y)\\}$ at quantile
$\alpha$ when $A=a_{0}$ is $\widehat{q}_{Y(a_{0})}^{{\bf
w}}(\alpha)=\inf\\{y:\alpha\leq\widehat{F}_{Y(a_{0})}^{{\bf w}}(y)\\}$.

The error of (5) at $A=a_{0}$ can be decomposed as

$\displaystyle\widehat{F}_{Y(a_{0})}^{{\bf w}}(y)-{F}_{Y(a_{0})}(y)={}$
$\displaystyle\int_{\mathcal{X}}F(y|{\bf
x},a_{0})\int_{\mathcal{A}}\mathrm{d}\left[F^{n}_{{\bf X},A,{\bf
w}}-F^{n}_{{\bf X}}F^{n}_{A}\right]({\bf x},a)$
$\displaystyle+R_{n}(a_{0})+R_{n,h}(a_{0})+R_{1n,h,{\bf
w}}(a_{0})+R_{2n,h,{\bf w}}(a_{0})$
$\displaystyle+\frac{1}{n}\sum_{i=1}^{n}e_{i}(y)w_{i}\widehat{f}^{-1}_{n,h}(a_{0}){K_{h}(A_{i}-a_{0})},$

where $e_{i}(y)=I(Y_{i}\leq y)-\mathbb{P}(Y\leq y|{\bf X}_{i},A_{i})$ where by
construction, $\mathbb{E}[e_{i}(y)]=0$, and $R_{n}(a_{0})$, $R_{n,h}(a_{0})$,
$R_{1n,h,{\bf w}}(a_{0})$, and $R_{2n,h,{\bf w}}(a_{0})$ behave similarly to
the second through fifth terms in the decomposition (2) with their consistency
for zero only depending on either the sample size and/or bandwidth $h$ and not
the weights. Thus, just as in estimation of the ADRF, the key source of
systematic bias in causal quantile function estimation can be mitigated by the
use of weights that make the covariates and treatment independent.

##  3 Measuring and Controlling Weighted Dependence with Energy Statistics

###  3.1 A criterion to evaluate the quality of a set of weights

Having established the relationship between the dependence of $A$ and ${\bf
X}$ and the error in a weighted nonparametric estimate of the ADRF, we now
construct a criterion that can assess how well a given set of weights induces
independence on the weighted scale, i.e., we aim to characterize and bound the
distance between $F^{n}_{{\bf X},A,{\bf w}}$ and $F^{n}_{{\bf X}}F^{n}_{A}$.
We do so by building on the ideas of distance covariance (Székely et al.,
2007), which is a measure of dependence between two random vectors of
arbitrary dimension. The population distance covariance is zero if and only if
the vectors are independent. Hence, a weighted distance covariance will be a
useful component for our measure. Let ${\bf w}=(w_{1},\dots,w_{n})$ and define
the weighted distance covariance to be

$\mathcal{V}^{2}_{n,{\bf w}}({\bf
X},A)=\frac{1}{n^{2}}\sum_{k,\ell=1}^{n}w_{k}w_{\ell}C_{k\ell}D_{k\ell},$ (6)

where

$\displaystyle c_{k\ell}$ $\displaystyle=\lVert{\bf X}_{k}-{\bf
X}_{\ell}\rVert_{2},\quad\mkern
1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu_{k\cdot}$
$\displaystyle=\frac{1}{n}\sum_{\ell=1}^{n}c_{k\ell},\quad\mkern
1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern
1.5mu_{\cdot\ell}=\frac{1}{n}\sum_{k=1}^{n}c_{k\ell},$ $\displaystyle\mkern
1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu_{\cdot\cdot}$
$\displaystyle=\frac{1}{n^{2}}\sum_{k,\ell=1}^{n}c_{k\ell},\quad C_{k\ell}$
$\displaystyle=c_{k\ell}-\mkern
1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu_{k\cdot}-\mkern
1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu_{\cdot\ell}+\mkern
1.5mu\overline{\mkern-1.5muc\mkern-1.5mu}\mkern 1.5mu_{\cdot\cdot},$

for $k,\ell=1,\dots,n$. Similarly define $d_{k\ell}=\lvert
A_{k}-A_{\ell}\rvert$, $\mkern 1.5mu\overline{\mkern-1.5mud\mkern-1.5mu}\mkern
1.5mu_{k\cdot}=\frac{1}{n}\sum_{\ell=1}^{n}d_{k\ell}$, $\mkern
1.5mu\overline{\mkern-1.5mud\mkern-1.5mu}\mkern
1.5mu_{\cdot\ell}=\frac{1}{n}\sum_{k=1}^{n}d_{k\ell}$, and
$D_{k\ell}=d_{k\ell}-\mkern 1.5mu\overline{\mkern-1.5mud\mkern-1.5mu}\mkern
1.5mu_{k\cdot}-\mkern 1.5mu\overline{\mkern-1.5mud\mkern-1.5mu}\mkern
1.5mu_{\cdot\ell}+\mkern 1.5mu\overline{\mkern-1.5mud\mkern-1.5mu}\mkern
1.5mu_{\cdot\cdot}$. The quantity (6) simplifies to the original distance
covariance when weights are all $1$. Since the original distance covariance is
always non-zero, (6) is also always non-zero if the weights are positive. We
now provide further insight and motivation of the form of
$\mathcal{V}^{2}_{n,{\bf w}}({\bf X},A)$ and its interpretation in terms of
weighted distributions.

Letting $i=\sqrt{-1}$, we define the (weighted) empirical characteristic
functions as $\varphi^{n}_{{\bf X},A,{\bf w}}({\bf t},s)=$
$\frac{1}{n}\sum_{i=1}^{n}w_{i}\exp\\{i\langle\mathbf{t},{{\bf
X}}_{i}\rangle+isA_{i}\\}$, $\varphi^{n}_{{\bf X},{\bf w}}({\bf t})=$
$\frac{1}{n}\sum_{i=1}^{n}w_{i}\exp\\{i\langle\mathbf{t},{{\bf
X}}_{i}\rangle\\}$, $\varphi^{n}_{A,{\bf
w}}(s)=\frac{1}{n}\sum_{i=1}^{n}w_{i}\exp\\{isA_{i}\\}$, and unweighted
empirical characteristic functions $\varphi^{n}_{{\bf X},A}({\bf t},s)$,
$\varphi^{n}_{{\bf X}}({\bf t})$ and $\varphi^{n}_{A}(s)$ are defined
accordingly.

######  Theorem 3.1.

Let ${\bf w}=(w_{1},\dots,w_{n})$ be a vector of weights such that
$\sum_{i=1}^{n}w_{i}=n$ and $w_{i}\geq 0$ for all $i=1,\dots,n$. Then
$\mathcal{V}^{2}_{n,{\bf w}}({\bf X},A)\geq 0$ and

$\displaystyle\mathcal{V}^{2}_{n,{\bf w}}({\bf X},A)={}$
$\displaystyle\int_{\mathbb{R}^{p+1}}|\varphi^{n}_{{\bf X},A,{\bf w}}({\bf
t},s)-\varphi^{n}_{{\bf X},{\bf w}}({\bf t})\varphi^{n}_{A,{\bf w}}(s)$
$\displaystyle\;\;\quad\quad+(\varphi^{n}_{{\bf X},{\bf w}}({\bf
t})-\varphi^{n}_{{\bf X}}({\bf t}))(\varphi^{n}_{A,{\bf
w}}(s)-\varphi^{n}_{A}(s))|^{2}\omega({\bf t},s)\mathrm{d}{\bf
t}\;\mathrm{d}s$ (7)

where $\omega({\bf t},s)=(c_{p}c_{1}\lVert{\bf
t}\rVert_{2}^{1+p}|s|^{2})^{-1}$ with
$c_{d}=\frac{\pi^{(1+d)/2}}{\Gamma((1+d)/2)}$, and $\Gamma(\cdot)$ is the
complete gamma function.

Based on (7), it is clear that if $F^{n}_{{\bf X},A,{\bf w}}=F^{n}_{{\bf
X}}F^{n}_{A}$, then $\mathcal{V}^{2}_{n,{\bf w}}({\bf X},A)=0$. However, the
converse is not necessarily true. Yet, if the weights preserve the marginal
distribution of treatment and covariates, i.e. $F^{n}_{{\bf X},{\bf
w}}=F^{n}_{{\bf X}}$ and $F^{n}_{A,{\bf w}}=F^{n}_{A}$, then
$\mathcal{V}^{2}_{n,{\bf w}}({\bf X},A)=0$ implies that $F^{n}_{{\bf X},A,{\bf
w}}=F^{n}_{{\bf X}}F^{n}_{A}$. In other words, if one can add additional terms
to (6) that also measure the distance between $F^{n}_{{\bf X},{\bf w}}$ and
$F^{n}_{{\bf X}}$ along with that between $F^{n}_{A,{\bf w}}$ and $F^{n}_{A}$,
then (6) can be utilized to construct a measure that determines the distance
between $F^{n}_{{\bf X},A,{\bf w}}$ and $F^{n}_{{\bf X}}F^{n}_{A}$, i.e., a
measure for the weighted dependence between ${\bf X}$ and $A$. In line with
our use of the distance covariance, we leverage the weighted energy distance
proposed in Huling and Mak (2020) to help construct such a measure.

Applied to our setting, the weighted energy distance between $F^{n}_{{\bf
X},{\bf w}}$ and $F^{n}_{\bf X}$ is

$\displaystyle\quad\mathcal{E}(F^{n}_{{\bf X},{\bf w}},F^{n}_{{\bf
X}})\equiv{}$
$\displaystyle\frac{2}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}w_{i}\lVert{{\bf
X}}_{i}-{{\bf
X}}_{j}\rVert_{2}-\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}w_{i}w_{j}\lVert{{\bf
X}}_{i}-{{\bf
X}}_{j}\rVert_{2}-\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\lVert{{\bf
X}}_{i}-{{\bf X}}_{j}\rVert_{2}.$

Due to Proposition 1 of Huling and Mak (2020), it can be shown that
$\mathcal{E}(F^{n}_{{\bf X},{\bf w}},F^{n}_{{\bf
X}})=\int_{\mathbb{R}^{p}}|\varphi^{n}_{\bf X}({\bf t})-\varphi^{n}_{{\bf
X},{\bf w}}({\bf t})|^{2}\omega(\mathbf{t})\mathrm{d}\mathbf{t}$, where
$\omega(\mathbf{t})=1/(C_{p}\lVert{\bf t}\rVert_{2}|^{1+p})$,
$C_{p}=\pi^{(1+p)/2}/\Gamma((1+p)/2)$ is a constant. The weighted energy
distance $\mathcal{E}(F^{n}_{A,{\bf w}},F^{n}_{A})$ between $F^{n}_{A,{\bf
w}}$ and $F^{n}_{A}$ can be similarly defined.

Our proposed measure of the level of independence induced by a set of weights
is defined as

${\cal D}({\bf w})=\mathcal{V}^{2}_{n,{\bf w}}({\bf
X},A)+\mathcal{E}(F^{n}_{{\bf X},{\bf w}},F^{n}_{{\bf
X}})+\mathcal{E}(F^{n}_{A,{\bf w}},F^{n}_{A}).$ (8)

The following result demonstrates that ${\cal D}({\bf w})$ indeed achieves its
stated goal.

######  Theorem 3.2.

Let ${\bf w}=(w_{1},\dots,w_{n})$ be a vector of weights with $w_{i}>0$ and
$\sum_{i=1}^{n}w_{i}=n$. Then ${\cal D}({\bf w})\geq 0$ with equality to zero
if and only if $\varphi^{n}_{{\bf X},A,{\bf w}}({\bf t},s)=\varphi^{n}_{{\bf
X}}({\bf t})\varphi^{n}_{A}(s)$, $\varphi^{n}_{{\bf X},{\bf w}}({\bf
t})=\varphi^{n}_{{\bf X}}({\bf t})$, and $\varphi^{n}_{A,{\bf
w}}(s)=\varphi^{n}_{A}(s)$ for all $({\bf t},s)\in\mathbb{R}^{p+1}$. Further,
$\int|\varphi^{n}_{{\bf X},A,{\bf w}}({\bf t},s)-\varphi^{n}_{{\bf X}}({\bf
t})\varphi^{n}_{A}(s)|^{2}\omega({\bf t},s)\mathrm{d}{\bf t}\mathrm{d}s\leq
3{\cal D}({\bf w})$.

Thus, smaller values of ${\cal D}({\bf w})$ indicate weaker dependence between
${\bf X}$ and $A$ after weighting and better preservation of the marginal
distributions and larger values implying the opposite. ${\cal D}({\bf w})=0$
implies the weights induce complete independence between $A$ and ${\bf X}$ and
that the marginal distributions of ${\bf X}$ and $A$ are exactly preserved.

We also have the following result, which shows how the proposed distance acts
as a bound on integration errors over a class of functions.

######  Lemma 3.3.

Let $\mathcal{H}$ be the native space induced by the radial kernel
$\Phi(\cdot,\cdot)=-\|\cdot\|_{2}\times-|\cdot|$ on
$\mathcal{X}\times\mathcal{A}$, and let $g(\cdot,\cdot)$ be any function where
$g(\cdot,\cdot)\in\mathcal{H}$. Then, for any weights ${\bf w}$ satisfying
$\sum_{i=1}^{n}w_{i}=n,w_{i}\geq 0$, we have

$\left[\int_{\mathcal{X}}\int_{\mathcal{A}}g({\bf
x},a)\mathrm{d}\left[F^{n}_{{\bf X},A,{\bf w}}-F^{n}_{\bf
X}F^{n}_{A}\right]({\bf x},a)\right]^{2}\leq C_{g}\mathcal{D}({\bf w}),$ (9)

where $C_{g}=3\|g\|^{2}_{\mathcal{H}}\geq 0$ is a constant depending on only
$g$.

For any $a\in\mathcal{A}$ if $\mu(\cdot,a)\in\mathcal{H}$, we can see that
${\cal D}({\bf w})$, modulo a constant, acts as a bound on the systematic bias
term (2) as long as $\mu(\cdot,a)$ is sufficiently smooth. Thus, if $\mu({\bf
x},a)$ is contained in $\mathcal{H}$, then we can expect weights with smaller
${\cal D}({\bf w})$ to lead to smaller systematic bias. The space
$\mathcal{H}$ is a reasonably broad class of functions as it contains the
Sobolev space of functions with square-integrable functions with
$r<\lceil(p+1)/2\rceil$-th differentials (Mak and Joseph, 2018; Huling and
Mak, 2020). Our goal for the next section is to define weights that are
optimal in terms of our criterion. The weights that minimize ${\cal D}({\bf
w})$ will result in mitigation of the dependence of $A$ and ${\bf X}$ induced
by nonrandom selection into treatment.

It is natural to wonder whether it would instead be more appropriate to simply
define a distance as $\int|\varphi^{n}_{{\bf X},A,{\bf w}}({\bf
t},s)-\varphi^{n}_{{\bf X}}({\bf t})\varphi^{n}_{A}(s)|^{2}\omega({\bf
t},s)\mathrm{d}{\bf t}\mathrm{d}s$ and construct a relationship between this
distance and Euclidean norms computable from data. However, we have found that
the resulting quantity can be empirically problematic and unreliable. Further,
we have found that optimization of such a quantity cannot be achieved reliably
by existing algorithms and is thus not suitable for the proactive construction
of weights. Our proposed quantity, while more complicated, does not exhibit
any of these issues in the sense that it reliably measures dependence, and, as
we will demonstrate, it is straightforward to optimize with existing quadratic
programming software.

###  3.2 A new proposal: optimal independence weights

We define the distance covariance optimal weights (DCOWs) to be

$\displaystyle{\bf w}^{d}_{n}\in{}$ $\displaystyle\operatorname*{argmin}_{{\bf
w}=(w_{1},\dots,w_{n})}{\cal D}({\bf w})\text{ such that
}\sum_{i=1}^{n}w_{i}=n,\text{ and }w_{i}\geq 0\text{ for }i=1,\dots,n.$ (10)

Due to Theorem 3.2, the distance covariance optimal weights ${\bf w}^{d}_{n}$
are designed to minimize dependence between ${\bf X}$ and $A$ on the weighted
scale while keeping the weighted marginal distributions of ${\bf X}$ and $A$
close to those of the unweighted data. The constraint $\sum_{i=1}^{n}w_{i}=n$
ensures that $F^{n}_{{\bf X},A,{{\bf w}^{d}_{n}}}$ is a valid distribution
function. Since ${\cal D}({\bf w})$ tracks with the dependence induced by a
set of weights, the DCOWs can be thought of as optimal independence weights.

Although, as we will show later, the DCOWs result in consistent weighted dose-
response estimators due to their explicit targeting of the source of
systematic error, without additional constraints they may not guarantee
optimal convergence rates. Instead, a small change to our criterion to allow
for penalization of the squares of the weights allows for us to control the
variability of the weights without sacrificing their bias-reduction property.
This additional penalty is akin to focusing more directly on mean squared
error instead of bias and can be interpreted as discouraging the effective
sample size (ESS) of the weights from being too small. In particular, the ESS
is typically approximated as $(\sum_{i}w_{i})^{2}/\sum_{i}w_{i}^{2}$ (Kish,
1965), which, due to our constraints, is $4n^{2}/\sum_{i}w_{i}^{2}$ and is
precisely the inverse of our proposed penalty. Further, combining a “balance”
criterion with a means of weight variability mitigation is in line with the
recommendations of Chattopadhyay et al. (2020).

We now define the penalized distance covariance optimal weights (PDCOWs) to be

$\displaystyle{\bf w}^{pd}_{n}\in{}$
$\displaystyle\operatorname*{argmin}_{{\bf w}=(w_{1},\dots,w_{n})}{\cal
D}({\bf w})+\lambda\frac{1}{n^{2}}\sum_{i=1}^{n}w_{i}^{2}$
$\displaystyle\text{ such that }\sum_{i=1}^{n}w_{i}=n,w_{i}\geq 0\text{ for
}0<\lambda<\infty,i=1,\dots,n.$ (11)

Here, the tuning parameter $\lambda$ is any positive constant and can be
chosen by the user to achieve a desired ESS. A lemma provided in the
Supplement similar to Lemma 3.3 shows that the penalized version of our
criterion acts as a bound on the term in the left hand side of (9) plus the
squares of the weights, which is more akin to a bound on the root mean squared
error than bias as in Lemma 3.3. Although additional penalty terms are
necessary for convergence rate guarantees, in practice we have found that only
very minimal or even no penalization at all works well because the unpenalized
weights, the DCOWs, tend to be quite stable. In all analyses described later,
we use only the DCOWs with no penalization of the weights. Alternatively, one
can simply choose $\lambda$ so that a desired effective sample size is
achieved.

Both the DCOWs and PDCOWs can be used in a wide variety of estimators for
various causal estimands involving continuous treatments, not just Nadaraya-
Watson-based estimators and not just estimators of the ADRF. The weights can
be used either in a simple weighted nonparametric estimator of the dose-
response function or to supplement any doubly-robust estimator of such. For
the former, our weights provide a completely nonparametric and empirically
robust estimation approach that requires only mild moment conditions on the
covariates and treatment for estimation consistency. For the latter, such an
estimator using our weights is guaranteed to be consistent regardless of the
correctness of the outcome model, while it still enjoys efficiency gains if
the outcome model is well-specified.

In a later section, we provide more formal statements on the consistency of
dose-response function estimators that use our proposed distance covariance
optimal weights.

###  3.3 Doubly robust estimation with optimal independence weights

Another class of estimators for causal ADRFs are doubly-robust estimators such
as in Kennedy et al. (2017) and Colangelo and Lee (2020), which combine sample
weights and an outcome model $\widehat{\mu}({\bf x},a_{0})$, ideally a
consistent estimator of ${\mu}({\bf x},a_{0})$. Though the term “doubly-
robust” is reflective of the property that if either the weights or outcome
model is correctly specified, the estimator will be consistent, a more
consequential property of doubly-robust estimators is that the error rates of
each model are multiplied. Loosely-speaking, this implies that reasonable
estimation of one model can alleviate the strain on the other. However, the
variance of such estimators tends to be more significantly impacted by
variable or poor estimation of the weights, indicating potentially an outsize
impact of poor estimation of the weights in practice relative to that of the
outcome model.

Although in the next section we reveal that the use of the DCOWs alone results
in the ideal convergence rate for the ADRF, our view is that DCOWs can be
enhanced by using doubly-robust estimators, or more provocatively, that the
use of DCOWs can significantly enhance doubly-robust estimators. The DCOWs
provide grounding and confidence to the analyst that the estimator will
converge at the right rate without the need for careful modeling, but the use
of a reasonably well-specified outcome model provides an opportunity to fine-
tune performance, resulting in an estimator that works well empirically.
Unlike in Colangelo and Lee (2020), we do not require the use of flexible
machine learning methods to guarantee consistency due to the independence-
inducing property of our weights.

Here, for simplicity of presentation, we focus on the following Nadaraya-
Watson-based doubly-robust estimator based on any estimator $\hat{\mu}({\bf
x},a_{0})$ of $\mu({\bf x},a_{0})$ as

$\mu_{NW}^{{\bf w},DR}(a_{0})=\frac{1}{n}\sum_{i=1}^{n}\hat{\mu}({\bf
X}_{i},a_{0})+\frac{\sum_{i=1}^{n}(Y_{i}-\hat{\mu}({\bf
X}_{i},a_{0}))w_{i}K_{h}(A_{i}-a_{0})}{\sum_{i=1}^{n}K_{h}(A_{i}-a_{0})}.$

In the Appendix we derive a decomposition of the error $\mu_{NW}^{{\bf
w},DR}(a_{0})-\mu(a_{0})$ and show that the systematic bias term not related
to smoothing is

$\int_{\mathcal{X}}\left\\{\mu({\bf x},a_{0})-\widehat{\mu}({\bf
x},a_{0})\right\\}\int_{\mathcal{A}}\mathrm{d}\left[F^{n}_{{\bf X},A,{\bf
w}}-F^{n}_{{\bf X}}F^{n}_{A}\right]({\bf x},a).$ (12)

Lemma 3.3 implies that (12) is less than or equal to $3\|\mu({\bf
x},a_{0})-\widehat{\mu}({\bf x},a_{0})\|_{\mathcal{H}}{\cal D}({\bf w})$
provided that $\mu(\cdot,a_{0})-\widehat{\mu}(\cdot,a_{0})\in\mathcal{H}$.
With the DCOWs, we provide a set of weights ${\bf w}^{d}$ that makes ${\cal
D}({\bf w})$ as small as possible, though it may not be exactly zero for a
finite sample.

Lemma 3.3 and our error decomposition for doubly-robust estimators illustrate
that a reasonable but imperfect model can help mitigate leftover bias due to
residual dependence, i.e. locations where $F^{n}_{{\bf X},A,{\bf
w}^{d}}-F^{n}_{{\bf X}}F^{n}_{A}$ differs significantly from zero. For a
doubly-robust estimator with a well-specified outcome model, the constant
$3\|g\|_{\mathcal{H}}$ in the bound of Lemma 3.3 can be smaller than for a
non-doubly robust estimator. As a result, using a doubly-robust estimator with
the weights ${\bf w}^{d}$ and a reasonably well-specified outcome model
results in a tighter bound on the systematic bias term. Further, the
additional terms in the error decomposition for $\mu_{NW}^{{\bf w},DR}(a_{0})$
that involve $\widehat{\mu}({\bf x},a_{0})$ are likely to be made smaller as
long as $\mu({\bf x},a_{0})\approx\widehat{\mu}({\bf x},a_{0})$. Given this,
it is clear that using a doubly-robust estimator with the PDCOWs is a well-
motivated strategy to improve performance for a given dataset. Although this
exploration focuses on a Nadaraya-Watson estimator, the same message applies
to the local linear doubly-robust estimator of Kennedy et al. (2017).

###  3.4 Asymptotic properties

We first show that the distance covariance optimal weights do indeed induce
independence asymptotically. Throughout this section we define $w({\bf
x},a)=f_{A}(a)/f_{A|{\bf X}}(a|{\bf X}={\bf x})$ to be the “true” normalized
density weights.

######  Theorem 3.4.

Let ${\bf w}^{d}_{n}$ be the distance covariance optimal weights defined in
(10). Then if $\mathbb{E}\lVert{\bf X}\rVert_{2}<\infty$ and $\mathbb{E}\lvert
A\rvert<\infty$ we have

$\lim_{n\to\infty}F^{n}_{{\bf X},A,{\bf w}^{d}}({\bf x},a)=F_{{\bf X}}({\bf
x})F_{A}(a)$ (13)

almost surely for every continuity point $({\bf x},a)\in\mathbb{R}^{p+1}$. If,
additionally $\mathbb{E}w^{2}({\bf X},A)<\infty$ holds, then

$\lim_{n\to\infty}F^{n}_{{\bf X},A,{\bf w}^{pd}}({\bf x},a)=F_{{\bf X}}({\bf
x})F_{A}(a)$ (14)

almost surely for every continuity point $({\bf x},a)\in\mathbb{R}^{p+1}$.

Theorem 3.4 is in some sense the most important property of the DCOWs, as it
demonstrates the feasibility of using these weights not just in dose-response
function estimation, but also in the estimation of any causal estimand that
requires independence.

We now show that for the particular task of estimating the ADRF, using the
DCOWs in a weighted Nadaraya-Watson estimator results in consistent estimation
of the dose-response function.

######  Theorem 3.5.

Assume that the kernel $K(\cdot)$ is symmetric, second order, i.e. it meets
the conditions that $\int uK(u)\mathrm{d}u=0$, $\int K(u)\mathrm{d}u=1$, and
$0<\int u^{2}K(u)\mathrm{d}u<\infty$, and is bounded differentiable. Further,
assume that the moment conditions required in Theorem 3.4 hold and that
$\mu({\bf x},a_{0})$ and is bounded and continuous on
$\mathcal{X}\times\mathcal{A}$ and has second order derivatives,
$f_{A}(a_{0})$ is bounded and has second order derivatives, $1/f_{A}(a_{0})$
is uniformly bounded. Then for ${\bf w}={\bf w}^{d}\text{ and }{\bf w}={\bf
w}^{pd}$

$\lim_{h\to 0,nh\to\infty}\mu_{NW}^{{\bf w}}(a_{0})=\mu(a_{0})$ (15)

in probability for all continuity points $a_{0}\in\mathcal{A}$.

Thus, both the DCOWs and PDCOWs result in consistent estimation of the causal
average dose-response function. Since $\mu_{NW}^{{\bf w}}(a_{0})$ can be
viewed as a special case of the double-robust estimator with a misspecified
outcome model, it can also be straightforwardly shown that $\mu_{NW}^{{\bf
w},DR}(a_{0})$ is also consistent even if $\widehat{\mu}({\bf x},a_{0})$ is
inconsistent for $\mu({\bf x},a_{0})$ as long as the outcome model does not
diverge.

The following result shows the rate of convergence of ${\cal D}({\bf w}^{d})$
and ${\cal D}({\bf w}^{pd})$ to zero and also shows that the squares of the
PDCOWs are controlled asymptotically under certain moment conditions.

######  Lemma 3.6.

Assume the moment conditions (A1) and (A2) listed in Appendix B hold, then
${\cal D}({\bf w}^{d})=O_{p}(1/n)$. If in addition we have
$\mathbb{E}w^{2}({\bf X},A)<\infty$, then ${\cal D}({\bf w}^{pd})=O_{p}(1/n)$
and $\frac{1}{n^{2}}\sum_{i=1}^{n}{{w_{i}^{pd}{}}^{2}}=O_{p}(1/n)$.

Lemma 3.6 implies the rate of convergence of $\int|\varphi^{n}_{{\bf X},A,{\bf
w}}({\bf t},s)-\varphi^{n}_{{\bf X}}({\bf
t})\varphi^{n}_{A}(s)|^{2}\omega({\bf t},s)\mathrm{d}{\bf t}\mathrm{d}s$ to
zero for both ${\bf w}^{d}$ and ${\bf w}^{pd}$, indicating they result in
independence of the treatment and covariates at a parametric rate. It also
implies that the variability of ${\bf w}^{pd}$ is not a rate-limiting factor
in nonparametric estimates of the ADRF.

The following shows the convergence rate of the penalized distance covariance
optimal weights under additional mild moment conditions on the covariates,
treatment, and $w({\bf X},A)$.

######  Theorem 3.7.

Assume the conditions required in Theorem 3.5 hold. Additionally assume the
moment conditions required in Lemma 3.6 hold and that
$\mathbb{E}[Y^{2}\>|\>{\bf X}={\bf x},A=a]<c$ for some $c$ uniformly over
${\bf x}\in\mathcal{X}$. Then

$\mu_{NW}^{{\bf w}^{pd}}(a_{0})-\mu(a_{0})=O_{p}(1/\sqrt{nh}+h^{2}).$ (16)

As long as the outcome model does not diverge, a result similar to Theorem 3.7
holds for the doubly-robust estimator $\mu_{NW}^{{\bf w},DR}(a_{0})$. This
convergence rate is the standard rate for unweighted Nadaraya-Watson
estimators of a univariate regression function; thus, the convergence rate of
the weighted estimator based on our proposed weights is unaffected by the
nonparametric nature of the estimation of ${\bf w}^{pd}$.

###  3.5 Computation

The optimization problems (10) and (11) can be formulated as quadratic
programming problems with linear constraints, making them straightforward to
implement with commercial and open-source solvers such as OSQP (Stellato et
al., 2020) or CCCP (Andersen et al., 2011). We will show how the (P)DCOW
optimization problem can be represented as a generic quadratic program

$\min_{{\bf w}}\>{\bf w}^{{\sf\scriptscriptstyle{T}}}{\bf P}{\bf w}+{\bf
q}^{{\sf\scriptscriptstyle{T}}}{\bf w}\ \text{ s.t. }\ {\bf
A}^{{\sf\scriptscriptstyle{T}}}{\bf w}\leq{\bf c}.$

The (P)DCOW optimization criterion takes ${\bf P}={\bf P}_{DCOW}+{\bf
P}_{EB(A)}+{\bf P}_{EB({\bf X})}+{\bf P}_{\lambda}$, ${\bf
q}^{{\sf\scriptscriptstyle{T}}}={\bf
q}^{{\sf\scriptscriptstyle{T}}}_{EB(A)}+{\bf
q}^{{\sf\scriptscriptstyle{T}}}_{EB({\bf X})}$, where ${\bf
P}_{DCOW}=\frac{1}{n^{2}}{\bf C}\circ{\bf D}$, ${\bf
P}_{EB(A)}=-\frac{1}{n^{2}}{\bf Q}_{A},{\bf P}_{EB({\bf
X})}=-\frac{1}{n^{2}}{\bf Q}_{X},{\bf P}_{\lambda}=\lambda\frac{1}{n^{2}}{\bf
I}_{n}$, ${\bf
q}^{{\sf\scriptscriptstyle{T}}}_{EB(A)}=\frac{2}{n^{2}}\mathbf{1}^{{\sf\scriptscriptstyle{T}}}{\bf
Q}_{A}$, and ${\bf q}^{{\sf\scriptscriptstyle{T}}}_{EB({\bf
X})}=\frac{2}{n^{2}}\mathbf{1}^{{\sf\scriptscriptstyle{T}}}{\bf Q}_{{\bf X}}$.
The terms ${\bf C}=(C_{k\ell})$ and ${\bf D}=(D_{k\ell})$ are as defined in
(6), $\circ$ is the elementwise product, and ${\bf Q}_{A}$ and ${\bf Q}_{{\bf
X}}$ are the Euclidean distance matrices of $A$ and ${\bf X}$, respectively.
From this formulation, we have $\mathcal{V}^{2}_{n,{\bf w}}({\bf X},A)={\bf
w}^{{\sf\scriptscriptstyle{T}}}{\bf P}_{DCOW}{\bf w}$,
$\mathcal{E}(F^{n}_{A,{\bf w}},F^{n}_{A})={\bf
w}^{{\sf\scriptscriptstyle{T}}}{\bf P}_{EB(A)}{\bf w}+{\bf
q}^{{\sf\scriptscriptstyle{T}}}_{EB(A)}{\bf w}+c_{A}$,
$\mathcal{E}(F^{n}_{{\bf X},{\bf w}},F^{n}_{{\bf X}})={\bf
w}^{{\sf\scriptscriptstyle{T}}}{\bf P}_{EB(X)}{\bf w}+{\bf
q}^{{\sf\scriptscriptstyle{T}}}_{EB(X)}{\bf w}+c_{X}$, and
$\lambda\frac{1}{n^{2}}\sum_{i=1}^{n}w_{i}^{2}={\bf
w}^{{\sf\scriptscriptstyle{T}}}{\bf P}_{\lambda}{\bf w}$, where $c_{A}$ and
$c_{X}$ are constants, so minimizing ${\bf w}^{{\sf\scriptscriptstyle{T}}}{\bf
P}{\bf w}+{\bf q}^{{\sf\scriptscriptstyle{T}}}{\bf w}$ amounts to minimizing
${\cal D}({\bf w})+\lambda\frac{1}{n^{2}}\sum_{i=1}^{n}w_{i}^{2}$, the
solution to which are the PDCOWs (or the DCOWs when $\lambda=0$). The
constraints ${\bf A}^{{\sf\scriptscriptstyle{T}}}$ and ${\bf c}$ can be
specified as ${\bf A}^{{\sf\scriptscriptstyle{T}}}=(-{\bf
I}_{n},\mathbf{1},-\mathbf{1})^{{\sf\scriptscriptstyle{T}}}$ and ${\bf
c}=(\mathbf{0}^{{\sf\scriptscriptstyle{T}}},n,-n)^{{\sf\scriptscriptstyle{T}}}$
so that ${\bf A}^{{\sf\scriptscriptstyle{T}}}{\bf w}\leq{\bf c}$ implies
$w_{i}\geq 0$ for all $i$ and $\sum_{i}w_{i}=n$. Although the optimization
problem has the form of a linearly constrained quadratic program, the problem
is not convex because ${\bf P}$ is not positive semi-definite. Despite this,
algorithms such as OSQP or CCCP can reliably provide high quality solutions to
the optimization problem.

###  3.6 Extensions and implementation considerations

If in practice additional emphasis on correlations of particular moments of
covariates and the treatment is of importance, the optimization criterion in
(10) can be straightforwardly modified to optimize ${\cal D}({\bf w})$ subject
to the constraints

$\sum_{i=1}^{n}w_{i}\left\\{m_{\bf X}({\bf X}_{i})-\overline{m}_{\bf X}({\bf
X}_{i})\right\\}\left\\{m_{A}(A_{i})-\overline{m}_{A}(A_{i})\right\\}^{{\sf\scriptscriptstyle{T}}}={\bf
0},$

where ${m}_{\bf X}$ is a vector of moments of ${\bf X}$, $\overline{m}_{\bf
X}$ is the average of the moments over the sample, ${m}_{A}$ is a vector of
moments of $A$, and similarly $\overline{m}_{A}$ is their sample average.
These constraints allow for the user to construct weights that minimize
dependence between ${\bf X}$ and $A$ while emphasizing decorrelation of
desired moments. If exactly zero correlation is infeasible, a constraint that
the marginal correlations are less than a threshold can be used instead,
similar to Zubizarreta (2015). Additional emphasis on marginal correlation may
improve finite sample performance if known moments are especially imbalanced
or highly related to outcomes. These exact or inexact constraints can be
straightforwardly incorporated into the quadratic programming problem
described above.

When the dimension of ${\bf X}$ is high, the relative emphasis of the energy
distance terms for $A$ and ${\bf X}$ in ${\cal D}({\bf w})$ can be imbalanced,
with less emphasis on the individual components of ${\bf X}$ relative to $A$.
To help mitigate this issue, one can instead use a version of ${\cal D}({\bf
w})$ that aims to re-balance the terms to be more comparable on a per-unit
basis like the following: ${\cal D}_{c}({\bf w})=\mathcal{V}^{2}_{n,{\bf
w}}({\bf X},A)+c_{{\bf X}}\mathcal{E}(F^{n}_{{\bf X},{\bf w}},F^{n}_{{\bf
X}})+c_{A}\mathcal{E}(F^{n}_{A,{\bf w}},F^{n}_{A})$, where $c_{{\bf
X}}+c_{A}=1$ and $c_{{\bf X}}/c_{A}=\sqrt{p}$. This modification has no
meaningful impact on the theoretical properties of ${\cal D}_{c}({\bf w})$, as
our proof of Theorem 3.2 applies to this modified distance.

##  4 Numerical experiments

###  4.1 Comparator methods

We use the following methods to estimate weights. We use a naïve method which
uses weights equal to identity (unweighted). We use stabilized GPS weights
computed four ways: a linear regression model for estimating
$\mathbb{E}(A|{\bf X})$ (i.e., the conditional mean of dose given covariates)
and a normal conditional density (“GPS normal”); a gamma regression model for
estimating $\mathbb{E}(A|{\bf X})$ and a gamma conditional density (“GPS
gamma”); a generalized boosted model for estimating $\mathbb{E}(A|{\bf X})$,
where the number of trees was chosen to minimize the weighted average absolute
correlation between the treatment and covariates as in Zhu et al. (2015), and
a normal conditional density (“GBM”); and a Bayesian Additive Regression Trees
model for estimating $\mathbb{E}(A|{\bf X})$ and a normal conditional density
(“BART”). For methods that estimate weights by directly inducing a lack of
correlation between moments of the treatment and covariates, we use the
covariate balancing generalized propensity score of Fong et al. (2018) and the
entropy balancing approach of Tübbicke (2020); Vegetabile et al. (2021). Among
this class of methods, we only consider these two as other approaches to
estimating weights that decorrelate pre-specified moments behave largely
similarly to each other (Tübbicke, 2020; Vegetabile et al., 2021). We use the
exactly-identified version of the covariate balancing generalized propensity
score (“CBPS”). We use versions of the entropy balancing approach that
decorrelate either all first order moments (“Entropy (1)”), all first order
moments and squared terms in continuous covariates (“Entropy (2)”), or all
first order moments, pairwise interactions, and squared terms in continuous
covariates (“Entropy (int)”). Any resulting weights greater than 500 when
standardized to sum to $n$ are truncated at 500. We use our proposed DCOWs
(“DCOW”) and the proposed DCOWS where we further induce exact decorrelation of
first order moments as discussed in Section 3.6 (“DCOW (dm)”), both using the
dimension adjustment described in Section 3.5. For each method, the weights
are used in a weighted local linear regression used to estimate the ADRF. For
the GPS (normal), GPS (gamma), DCOW, and DCOW (dm) methods, we also use the
doubly-robust estimator of Kennedy et al. (2017) with a an outcome model that
is linear in the covariates with additive first order terms; the methods are
labeled as “GPS (normal,DR)”, “GPS (gamma,DR)”, “DCOW (DR)”, and “DCOW
(dm,DR)”, respectively.

###  4.2 Simulation using National Medical Expenditure Survey Data

The National Medical Expenditure Survey (NMES) relates medical expenditures
with degree of smoking among U.S. citizens (Johnson et al., 2003). The NMES
dataset contains information on 9708 individuals. The outcome variable is the
total medical expenditures in dollars and the treatment $A\in[0.05,216]$ is
the amount of smoking in pack years. We limit the data to those with $A\leq
80$ (i.e. $\mathcal{A}=[0.05,80]$), leaving 9368 units. Two of the covariates
are continuous ($X_{1}$: age started to smoke and $X_{2}$: age last smoked)
and the remaining 7 are discrete ($X_{3k}$: gender (male, female); $X_{4k}$:
race (African American, other); $X_{5k}$: seat belt use, (low, medium, high);
$X_{6k}$: education ($<$ high school, high school, $1-3$ years of college,
$4+$ years of college); $X_{6k}$: marital status (never married,
widowed/divorced, married); $X_{7k}$: census region (Northeast, Midwest,
South, West); $X_{8k}$: poverty status (poor, near poor, low income, middle
income, high income), with anywhere from $\ell_{k}=$ 2 to 5 levels. The
overall covariate vector thus has dimension 18. In our simulation, we leave
the treatment level (pack years) and covariates intact and simulate outcomes
for each unit from the following model.

$\displaystyle Y={}$
$\displaystyle\gamma_{1}b_{1}X_{1}+\gamma_{2}b_{2}\widetilde{X}_{1}^{2}+\gamma_{3}b_{3}\widetilde{X}_{1}^{3}+\gamma_{4}b_{4}X_{2}+\gamma_{5}b_{5}\widetilde{X}_{2}^{2}+\gamma_{6}b_{6}\widetilde{X}_{2}^{3}$
$\displaystyle+\sum_{j=3}^{8}\sum_{k=1}^{\ell_{j}}\left\\{\alpha_{jk}b_{1jk}+\eta_{1jk}b_{2jk}X_{1}+\eta_{2jk}b_{3jk}\widetilde{X}_{1}^{2}+\eta_{3jk}b_{4jk}X_{2}+\eta_{4jk}b_{5jk}\widetilde{X}_{2}^{2}\right\\}I(X_{j}=k)$
$\displaystyle-\overline{m}({\bf X})+f(A)(1+\delta({\bf X})))+\varepsilon,$

where $\widetilde{X}_{j}=X_{j}-\overline{X}_{j}$,
$\gamma_{1},\gamma_{3},\eta_{1jk},\eta_{3jk}\sim$ Unif$(-0.5,0.5)$,
$\gamma_{2},\gamma_{4},\eta_{2jk},\eta_{4jk}\sim$ Unif$(-0.1,0.1)$,
$\gamma_{3},\gamma_{6}\sim$ Unif$(-0.01,0.01)$, $\alpha_{jk}\sim$
Unif$(-10,10)$, $b_{j},b_{ijk}\sim$ Bernoulli$(0.5)$, $\overline{m}({\bf X})$
is the mean of the main effect terms preceding it, and the treatment effect
curve $f(A)=A/4+\frac{2}{(A/100+1/2)^{3}}-(A-40)^{2}/100$. For the constant
treatment effect setting, $\delta({\bf X})=0$, and for the heterogeneous
treatment effect setting,

$\delta({\bf
X})=\gamma_{I1}\widetilde{X}_{1}+\gamma_{I2}\widetilde{X}_{2}+\sum_{j=3}^{8}\sum_{k=1}^{\ell_{j}}\left\\{\eta_{I1jk}b_{I1jk}\widetilde{X}_{1}+\eta_{I2jk}b_{I2jk}\widetilde{X}_{2}\right\\}\widetilde{X}_{j,k}$

where $\widetilde{X}_{j,k}=I(X_{j}=k)-\overline{I(X_{j}=k)}$,
$\eta_{I1jk},\eta_{I2jk}\sim$ Unif$(-0.5,0.5)$, and $b_{I1jk},b_{I3jk}\sim$
Bernoulli$(0.5)$. We generate 100 different draws of the coefficients ($\eta$,
$\gamma$, and $b$ terms) in the outcome model above, allowing for the
simulation study to explore a wide variety of outcome models. For each of
these 100 outcome model draws, we replicate the simulation experiment 1000
times. For each replication, a random subsample of size $n<9368$ is drawn
without replacement from the 9368 units, and outcomes for these units are
generated from the outcome model. The simulation process is repeated for each
of the sample sizes $n=100,200,400,800,1600,3200$.

We also consider a simulation setting with the heterogeneous treatment effect
where 50 additional variables are added to the covariate vector so that the
overall dimension is 68. The 50 additional variables are generated so that
they are correlated with both the response and treatment while preserving the
original values of the response and treatment. To generate the additional
variables, we first draw random variables from a normal distribution with mean
$=(\overline{X}_{1}+\overline{X}_{2})/2$ and variance
$=(\text{var}(X_{1})+\text{var}(X_{2}))/2$, where $X_{1}$ is the age started
to smoke and $X_{2}$ is the age last smoked. We then use a Cholesky
decomposition to induce an $AR(1)$ correlation structure with correlation
parameter $\rho=-1/3$, where the correlation structure is with respect to the
expanded vector $(Y,A,X_{nv1},\dots,X_{nv50})$, where $X_{nvj}$ are the added
variables and the $AR(1)$ correlation structure is modified in such a way that
$A$ and $Y$ remain unchanged but the added variables have strong correlations
jointly with $A$ and $Y$ and thus impact the finite sample error in estimating
the ADRF if the added variables are not properly adjusted for. This variance
structure is used so that the variation of the noise variables is comparable
to that of the continuous covariates in the NMES data.

All methods are evaluated with a measure of the mean absolute bias (MAB) and
integrated root mean squared error (IRMSE), both of which are used for
evaluation of estimates of the ADRF in Kennedy et al. (2017). MAB and IRMSE
are defined as

$\displaystyle\text{MAB}={}$
$\displaystyle\int_{\mathcal{A}}\left\lvert\frac{1}{S}\sum_{s=1}^{S}\widehat{\mu}_{s}(a)-\mu(a)\right\rvert\widehat{f}_{A}(a)\mathrm{d}a\text{
and}$ $\displaystyle\text{IRMSE}={}$
$\displaystyle\int_{\mathcal{A}}\left[\frac{1}{S}\sum_{s=1}^{S}\left\\{\widehat{\mu}_{s}(a)-\mu(a)\right\\}^{2}\right]^{1/2}\widehat{f}_{A}(a)\mathrm{d}a,$

where $\widehat{f}_{A}(a)$ is a kernel density estimate of the marginal
density of the treatment variable, $s$ indexes the simulation replications,
and $\mathcal{A}$ is a trimmed version of the support of the treatment
variable that excludes pack years greater than 80. We calculate the MAB and
IRMSE statistics for each of the 100 different outcome model settings and then
compute their averages over the 100 settings.

The MAB and IRMSE results for the constant treatment effect setting are
displayed in Table 1 and the results for the heterogeneous treatment effect
setting are displayed in Table 2. Results for any method are not shown if no
numerical solution is found in more than 75% of the replications. The standard
DCOWs performed the best in terms of both MAB and IRMSE across all sample
sizes among all non-doubly-robust estimators, only being outperformed by the
doubly-robust estimator that uses DCOWs as weights. DCOWs with exact first
order moment decorrelation performed similarly to, but slightly worse than
standard DCOWs, for both the non-doubly-robust estimator and the doubly-robust
estimator, though the performance was much worse for $n=100$, as the exact
constraints may have been too stringent for the sample size. The estimators
using standard GPS weights, both doubly-robust and non-doubly-robust,
performed poorly in terms of both MAB and IRMSE for small to moderate sample
sizes, though the doubly-robust estimator with gamma regression-based weights
performed well in terms of MAB and IRMSE for larger sample sizes. The machine
learning approaches to GPS estimation (GBM and BART) performed relatively
poorly in terms of MAB and IRMSE for small sample sizes, with BART performing
better than GBM for small sample sizes, but with similar but slightly better
performance for GBM for large sample sizes. Among the GPS-based methods, GBM
and BART generally yielded the lowest MAB and IRMSE across the sample size
settings, except for the doubly-robust gamma GPS estimator. Among the moment
balancing approaches, entropy balancing with decorrelation induced only for
first and second order moments and not interactions performed the best in
terms of MAB and IRMSE, though for small sample sizes, entropy balancing
frequently failed to arrive at a solution. For the heterogeneous setting,
entropy balancing with decorrelation induced only for first order moments
performed the best, likely due to numerical instability of adding more moment
constraints. The CBPS estimator yielded worse performance than did entropy
balancing overall, though it did not face the convergence problems of entropy
balancing for small sample sizes. Entropy balancing with higher order moment
constraints failed to arrive at a solution far more frequently and performed
poorly even for large sample sizes, likely due to the strictness of the exact
moment constraints.

The results for the heterogeneous effect setting with 50 additional noise
variables are displayed in Table 3. Again, DCOWs performed the best among all
methods in terms of both MAB and IRMSE for all sample sizes except $n=100$,
where the doubly-robust gamma GPS estimate was slightly better in terms of
IRMSE.

Additional simulation results under the same setting of Vegetabile et al.
(2021) can be found in the supplementary file wherein our weights also yielded
the best performance across all sample size settings. The main difference was
that DCOWs without doubly-robust estimators performed slightly better than
DCOWs with doubly-robust estimators, likely due to substantial
misspecification of the outcome model.

$n=100$ $n=200$ $n=400$ $n=800$ $n=1600$ $n=3200$ Method MAB IRMSE MAB IRMSE
MAB IRMSE MAB IRMSE MAB IRMSE MAB IRMSE Unweighted 11.461 17.488 11.283 14.753
11.252 13.182 11.246 12.297 11.254 11.799 11.237 11.491 GPS (normal) 10.145
27.180 15.033 29.187 18.197 28.056 20.786 27.650 21.832 25.703 22.611 24.432
GPS (gamma) 9.340 23.994 9.892 19.269 9.993 15.831 10.112 13.978 10.084 12.693
10.006 11.162 GPS (normal,DR) 8.315 30.551 11.081 33.342 13.216 30.293 15.320
27.636 15.324 21.199 15.769 18.312 GPS (gamma,DR) 7.045 26.437 6.866 16.463
6.905 12.044 6.920 10.168 6.768 8.985 6.745 7.612 CBPS 9.871 25.398 12.463
23.375 14.055 22.016 14.064 18.802 14.970 17.546 15.712 16.847 GBM 8.190
23.428 7.551 17.298 7.105 12.498 7.379 10.432 7.524 9.225 7.746 8.589 BART
6.550 17.206 7.275 14.423 8.023 12.326 8.638 11.292 9.033 10.659 9.446 10.366
Entropy (1) —— —— 9.409 17.880 9.166 12.946 9.161 11.321 9.023 10.258 8.952
9.537 Entropy (2) —— —— 8.909 25.080 9.343 15.520 8.927 12.360 8.503 10.456
8.231 9.245 Entropy (2,int) —— —— —— —— —— —— —— —— 8.998 253.786 10.974
16.780 DCOW 4.684 12.204 3.866 8.383 3.495 6.245 3.252 4.947 2.992 4.075 2.750
3.416 DCOW (dm) 4.522 16.465 3.907 8.988 3.544 6.364 3.276 5.018 3.001 4.147
2.775 3.497 DCOW (DR) 3.904 9.284 3.335 6.455 2.663 4.586 2.196 3.388 1.919
2.661 1.753 2.189 DCOW (dm,DR) 3.905 11.567 3.459 6.919 2.828 4.806 2.356
3.568 2.057 2.817 1.873 2.320

Table 1: Mean absolute bias (MAB) and integrated root mean squared error
(IRMSE) for the constant treatment effect setting.

$n=100$ $n=200$ $n=400$ $n=800$ $n=1600$ $n=3200$ Method MAB IRMSE MAB IRMSE
MAB IRMSE MAB IRMSE MAB IRMSE MAB IRMSE Unweighted 20.075 28.198 19.865 24.382
19.863 22.361 19.826 21.165 19.860 20.552 19.838 20.151 GPS (normal) 16.490
36.634 20.904 36.562 23.999 34.649 26.269 33.482 27.298 31.335 27.982 29.873
GPS (gamma) 15.597 33.153 15.491 26.710 15.263 22.321 14.942 19.740 14.698
18.013 14.414 16.031 GPS (normal,DR) 12.062 41.168 15.197 40.612 17.611 36.146
19.466 32.248 19.562 25.719 19.958 22.685 GPS (gamma,DR) 12.021 38.538 11.194
25.270 10.856 18.518 10.553 15.730 10.341 14.197 10.235 11.970 CBPS 17.990
36.480 20.534 33.170 22.074 31.068 20.669 26.016 20.196 23.072 20.213 21.503
GBM 15.174 34.035 13.646 25.480 12.849 19.368 12.832 16.467 12.908 14.925
13.089 14.087 BART 13.372 26.339 13.141 21.772 13.753 18.962 14.286 17.460
14.712 16.654 15.118 16.245 Entropy (1) —— —— 12.066 24.966 11.030 16.764
10.545 13.917 10.229 12.200 10.061 11.053 Entropy (2) —— —— 15.737 36.405
15.872 24.570 14.401 19.456 13.181 16.191 12.242 13.928 Entropy (2,int) —— ——
—— —— —— —— —— —— 14.777 363.289 16.761 24.311 DCOW 6.718 18.517 6.083 13.254
5.818 10.552 5.569 8.852 5.400 7.703 5.329 6.830 DCOW (dm) 7.999 27.130 6.366
14.836 5.943 11.074 5.699 9.226 5.536 8.038 5.500 7.129 DCOW (DR) 5.778 16.726
5.220 12.206 4.891 9.780 4.735 8.311 4.829 7.422 5.055 6.753 DCOW (dm,DR)
6.447 21.900 5.569 13.414 5.158 10.321 5.007 8.760 5.088 7.823 5.276 7.074

Table 2: Mean absolute bias (MAB) and integrated root mean squared error
(IRMSE) for the heterogeneous treatment effect setting.

$n=100$ $n=200$ $n=400$ $n=800$ $n=1600$ $n=3200$ Method MAB IRMSE MAB IRMSE
MAB IRMSE MAB IRMSE MAB IRMSE MAB IRMSE Unweighted 20.075 28.198 19.865 24.382
19.863 22.361 19.826 21.165 19.860 20.552 19.838 20.151 GPS (normal) 19.733
29.907 16.453 31.114 22.881 39.654 30.521 46.508 37.341 50.345 42.592 51.453
GPS (gamma) 16.508 43.818 18.124 35.312 18.156 28.509 18.197 26.357 18.704
25.993 19.429 24.759 GPS (normal,DR) 13.993 20.345 14.330 21.161 17.472 33.781
22.528 42.084 26.955 41.828 30.583 39.910 GPS (gamma,DR) 13.655 74.504 15.441
46.099 15.694 26.679 15.769 22.863 16.455 22.531 17.155 20.910 CBPS 17.821
32.365 17.366 29.713 20.020 30.437 23.032 30.789 25.835 31.782 27.623 32.080
GBM 16.372 30.524 14.500 24.168 14.961 21.193 15.211 19.662 15.348 18.582
15.277 17.573 BART 14.451 25.927 13.831 23.202 15.749 22.567 17.324 22.518
17.971 21.877 17.587 20.575 Entropy (1) —— —— —— —— —— —— 14.089 18.085 13.380
15.623 13.123 14.327 Entropy (2) —— —— —— —— —— —— —— —— 18.963 28.307 17.360
19.206 Entropy (2,int) —— —— —— —— —— —— —— —— —— —— —— —— DCOW 12.589 24.710
11.362 18.520 11.310 15.067 11.435 13.282 11.275 12.169 10.825 11.239 DCOW
(dm) 20.054 28.215 14.120 31.635 13.428 20.513 12.623 15.718 11.289 12.633
10.117 10.691 DCOW (DR) 13.481 21.034 12.426 16.655 11.545 14.065 10.784
12.206 10.201 10.966 9.495 9.875 DCOW (dm,DR) 14.084 20.742 13.024 20.867
11.838 16.022 11.035 13.202 10.377 11.400 9.650 10.106

Table 3: Mean absolute bias (MAB) and integrated root mean squared error
(IRMSE) for the heterogeneous treatment effect setting with 50 noise variables
added.

##  5 Analysis of mechanical power data

We use the Medical Information Mart for Intensive Care III (MIMIC-III)
database Johnson et al. (2016) to study the impact of a large degree of
mechanical power of ventilation on mortality among critically ill patients in
an ICU using electronic health record (EHR) data. Our study and the
construction of the cohort from the MIMIC-III database is based the original
study of Neto et al. (2018) and the code provided by the authors located at
https://github.com/alistairewj/mechanical-power. As there are widely-used
formal guidelines that influence the management of ventilation among patients
with respiratory distress (Siegel et al., 2019; Papazian et al., 2019), there
are many observable factors in EHR data that are highly related to the
mechanical power of ventilation. Many of these are also closely related to
patient mortality. The guidelines involve consideration both of factors
individually and also of many interactions of these factors. In addition, as
the exposure, the mechanical power of ventilation, is itself a complex summary
of multiple manipulable elements of a ventilator, the guidelines are not
directly related to the exposure, but have a strong indirect effect on the
power of ventilation. Thus, the dependence between observable factors and the
exposure level is highly complex; methods to control for confounding in this
setting must be able to handle such complexity.

Patients included in the study were at least 16 years of age and received
invasive mechanical ventilation for at least 48 hours. The study contains 5014
patients, and the treatment variable of interest is the amount of energy
generated by the mechanical ventilator measured by the mean of the largest and
smallest mechanical power of ventilation in Joules per minute in the second 24
hour period in the ICU, as in Neto et al. (2018). The outcome is an indicator
of in-hospital mortality. As in the original analysis of Neto et al. (2018),
we limit our analysis to patients receiving less than or equal to 50 Joules
per minute, resulting in a sample size of 4933. 73 pre-treatment covariates,
some of which were discrete, were included to be balanced, resulting in a
total dimension of 97 of the vector ${\bf X}$ of potential confounders. All
methods considered in the simulation section were then applied to construct
weights aiming to control for dependence between the 97 covariates and
mechanical power. Both GPS approaches, GBM, and BART resulted in several
extraordinarily large weights; these weights were truncated to mitigate
extreme variation. We attempted to construct entropy balancing weights that
either exactly or approximately (within a correlation tolerance of $0.1$)
balance pairwise interactions, but no numerical solution was found.

For all methods, balance statistics, including our developed criterion (8) and
weighted correlations between first order moments and pairwise interactions of
covariates and mechanical power, are summarized in Table 4. In terms of
weighted marginal correlations, the version of our DCOWs that induces first
order marginal correlations between covariates and mechanical power to be zero
performs the best, with the standard DCOWs and entropy balancing weights a
close second. In terms of our proposed criterion (8), by definition the DCOWs
yield the smallest value. However it is notable that only a small price is
paid in terms of both (8) and effective sample size (ESS) in order to
additionally exactly decorrelate marginal covariate moments and treatment.
Among methods that do not target independence, entropy balancing has the
smallest value of (8), indicating that decorrelating first order moments in
this particular dataset mitigates a vast majority of the dependence between
covariates and mechanical power of ventilation. We also note that exactly
decorrelating second order moments using entropy balancing results in
additionally instability and thus overall worse mitigation of dependence.
However, exact decorrelation of marginal results in much larger standard
errors of the resulting ADRF, as evidenced in Figure 2. Interestingly, both in
terms of marginal weighted correlations and in terms of our proposed
independence metric, the flexible machine learning approaches (GBM and BART)
perform significantly worse than a simple normal model for the conditional
density.

| Unweighted |  | GPS
---
(normal)
| GPS
---
(gamma)
CBPS | GBM | BART |  | Entropy
---
(1)
| Entropy
---
(2)
DCOW |  | DCOW
---
(dm)
(8) | 21.0406 | 2.5735 | 4.2687 | 4.7127 | 12.9924 | 6.1664 | 1.0927 | 4.2621 | 0.2201 | 0.2366
${\cal D}({\bf w})$ | 21.0406 | 2.3492 | 4.0143 | 4.3284 | 12.9604 | 5.7097 | 1.0052 | 4.0832 | 0.1436 | 0.1512
$\mathcal{E}(F^{n}_{A,{\bf w}},F^{n}_{A})$ | 0.0000 | 0.0266 | 0.0031 | 0.0099 | 0.0030 | 0.0426 | 0.0296 | 0.0006 | 0.0508 | 0.0588
$\mathcal{E}(F^{n}_{{\bf X},{\bf w}},F^{n}_{{\bf X}})$ | 0.0000 | 0.1976 | 0.2513 | 0.3743 | 0.0290 | 0.4141 | 0.0579 | 0.1783 | 0.0258 | 0.0265
ESS | 4933 | 2048 | 1832 | 1017 | 3892 | 782 | 2185 | 478 | 2232 | 2057
mean($|$Corr$|$) | 0.0697 | 0.0142 | 0.0365 | 0.0246 | 0.0577 | 0.0452 | 0.0067 | 0.0208 | 0.0067 | 0.0029
sd($|$Corr$|$) | 0.0821 | 0.0151 | 0.0363 | 0.0249 | 0.0609 | 0.0434 | 0.0075 | 0.0219 | 0.0064 | 0.0031
median($|$Corr$|$) | 0.0405 | 0.0096 | 0.0260 | 0.0178 | 0.0361 | 0.0352 | 0.0046 | 0.0159 | 0.0051 | 0.0020
95-pctl($|$Corr$|$) | 0.2482 | 0.0439 | 0.1048 | 0.0633 | 0.2000 | 0.1177 | 0.0200 | 0.0547 | 0.0184 | 0.0085
max($|$Corr$|$) | 0.7877 | 0.2003 | 0.8536 | 0.4702 | 0.5682 | 0.9388 | 0.1225 | 0.5066 | 0.0839 | 0.0451

Table 4: Summary statistics (mean, standard deviation, median, 95th
percentile, and maximum) of the absolute weighted correlations of the first
five powers of mechanical power of ventilation and all marginal moments of
covariates, pairwise interactions of covariates, and up to 5th order
polynomials of covariates. Summaries are over $8284\times 5$ weighted
correlations of covariate moments and mechanical power moments.

For a concrete example of how DCOWs mitigate dependence, we illustrate the
unadjusted marginal dependence of the Pa$\text{O}_{2}$/Fi$\text{O}_{2}$ ratio
and the mechanical power of ventilation. Pa$\text{O}_{2}$/Fi$\text{O}_{2}$
ratio is a a pre-treatment covariate that characterizes acute hypoxemia,
defines presence and severity of acute respiratory distress syndrome (ADRS),
and plays a critical role in ventilation guidelines (Papazian et al., 2019).
The Pa$\text{O}_{2}$/Fi$\text{O}_{2}$ ratio is also strongly associated with
mortality and is thus a critical confounder in this study. We display weighted
marginal dependence of the Pa$\text{O}_{2}$/Fi$\text{O}_{2}$ ratio and the
mechanical power of ventilation with our proposed weights and with the
entropy-balancing weights in Figure 1. Compared with methods that aim to
exactly decorrelate specified moments, our proposed weights are able to handle
nonlinear dependence between covariates and treatment even in datasets with
moderately high dimension. While entropy balancing accounts for a significant
proportion of dependence, there remains residual nonlinear dependence.

Figure 1:  Shown are plots of the relationship between the minimum
Pa$\text{O}_{2}$/Fi$\text{O}_{2}$ ratio on day 1 in the ICU (on a log base 10
scale) and the treatment, including an unadjusted plot (left) and plots
adjusted by DCOWs and entropy balancing weights (right two plots). In the
adjusted plots, the transparency of each point is proportional to its assigned
weight, with lighter points indicating less weight. The blue line is a
weighted nonparametric regression of the treatment on
Pa$\text{O}_{2}$/Fi$\text{O}_{2}$ (on the log base 10 scale) and the red line
is a weighted linear regression.

For each method, we construct pointwise confidence bands of the ADRF using a
nonparametric bootstrap similar to Wang and Wahba (1995). Weighted local
linear regression estimates of the ADRF of mechanical power of ventilation on
in hospital mortality and 95% confidence bands are displayed in Figure 2. The
entropy balancing approach with higher order moments was not included due to
the small effective sample size and that large number of bootstrap
replications for which no numerical solution was found. While both the entropy
balancing weights and the GPS estimated with a normal model exhibit reasonable
control over dependence between covariates and treatment, they result in
unacceptably high variance in the estimate of the ADRF.

Figure 2:  Shown are weighted local linear regression estimates of the ADRF of
mechanical power on the probability of in hospital mortality and pointwise 95%
confidence bands estimated via a nonparametric bootstrap. 38 replications of
the bootstrap resulted in no numerical solution for entropy balancing weights
and are left out of the confidence interval calculation.

##  6 Discussion

In this paper we provide a detailed inspection of the role of weights in
weighted nonparametric estimates of causal quantities involving continuous-
valued exposures from observational data. This inspection shows clearly that
the key source of bias depends on the degree to which the weights induce
independence between the treatment and confounders. We then provide a measure
that characterizes how well a set of weights mitigates the dependence between
the treatment and confounders. This measure does not require any tuning
parameters, making it straightforward to deploy in practice. Given some light
smoothness conditions on the outcome data-generating model, this measure acts
as an upper bound on the key source of systematic bias in a weighted
nonparametric estimate of the ADRF. Our proposed weights, the DCOWS, which
minimize our measure of dependence, provide an empirically robust means for
estimating weights, as they directly target independence between the treatment
and confounders. These weights are a natural complement to doubly-robust
estimators, as they provide an anchor for the doubly-robust estimator: since
they are guaranteed to be consistent on their own, the consistency of the
doubly-robust estimator does not critically depend on the correctness of the
outcome model. Thus, the outcome model can be safely used as a tool by which
to reduce variability in the estimate of the ADRF.

Appendix

##  Appendix A Properties of stabilized generalized propensity score weights

The stabilized GPS weights $f_{A}(A)/f_{A|{\bf X}}(A\>|\>{\bf X})$ have
several key properties. Namely, when weighting by these weights they 1) result
in independence of $A$ and ${\bf X}$, 2) preserve the marginal distributions
of ${\bf X}$ and $A$, and 3) have mean 1. For the following, we let $g({\bf
X})$ be any function of ${\bf X}$ only and let $h(A)$ be any function of $A$
only. It can be straightforwardly shown that

$\displaystyle\mathbb{E}_{{\bf X},A}\left(g({\bf
X})h(A)\frac{f_{A}(A)}{f_{A|{\bf X}}(A\>|\>{\bf X})}\right)$
$\displaystyle=\mathbb{E}_{\bf X}(g({\bf X}))\mathbb{E}_{A}(h(A)),$
$\displaystyle\mathbb{E}_{{\bf X},A}\left(g({\bf X})\frac{f_{A}(A)}{f_{A|{\bf
X}}(A\>|\>{\bf X})}\right)={}\mathbb{E}_{{\bf X}}\left(g({\bf X})\right),$
$\displaystyle\text{ and }\mathbb{E}_{{\bf
X},A}\left(h(A)\frac{f_{A}(A)}{f_{A|{\bf X}}(A\>|\>{\bf
X})}\right)={}\mathbb{E}_{A}\left(h(A)\right).$

The first equality shows that when weighting by $\frac{f_{A}(A)}{f_{A|{\bf
X}}(A\>|\>{\bf X})}$, $A$ and ${\bf X}$ are independent, as the above holds
for any $g$ and $h$. Further, the second two equalities show that, on the
weighted scale, the marginal distributions of ${\bf X}$ and $A$ are preserved.
Finally, the weights have mean 1, i.e. $\mathbb{E}_{{\bf
X},A}\left(\frac{f_{A}(A)}{f_{A|{\bf X}}(A\>|\>{\bf X})}\right)=1$.

##  Appendix B Regularity conditions

For $({\bf X}_{1},A_{1}),\dots,({\bf
X}_{6},A_{6})\mathbin{\overset{i.i.d.}{\kern 0.0pt\sim}}F_{{\bf X},A}$, define

$k(({\bf X}_{1},A_{1}),\dots,({\bf X}_{6},A_{6}))=w({\bf X}_{1},A_{1})g_{\bf
X}({\bf X}_{1},{\bf X}_{2},{\bf X}_{3},{\bf
X}_{4})g_{A}(A_{1},A_{2},A_{5},A_{6})w({\bf X}_{2},A_{2})$

with $w({\bf x},a)=\frac{f_{A}(a)}{f_{A|{\bf X}}(a|{\bf X}={\bf x})}$,

$\displaystyle g_{\bf X}({\bf X}_{1},{\bf X}_{2},{\bf X}_{3},{\bf X}_{4})={}$
$\displaystyle\|{\bf X}_{1}-{\bf X}_{2}\|_{2}-\|{\bf X}_{1}-{\bf
X}_{3}\|_{2}-\|{\bf X}_{2}-{\bf X}_{4}\|_{2}+\|{\bf X}_{3}-{\bf
X}_{4}\|_{2},\text{ and}$ $\displaystyle g_{A}(A_{1},A_{2},A_{3},A_{4})={}$
$\displaystyle|A_{1}-A_{2}|-|A_{1}-A_{3}|-|A_{2}-A_{4}|+|A_{3}-A_{4}|.$

Further define

$\displaystyle k_{\bf X}(({\bf X}_{1},A_{1}),\dots,({\bf X}_{4},A_{4}))={}$
$\displaystyle w({\bf X}_{1},A_{1})\|{\bf X}_{1}-{\bf X}_{3}\|_{2}+w({\bf
X}_{2},A_{2})\|{\bf X}_{2}-{\bf X}_{4}\|_{2}$ $\displaystyle-w({\bf
X}_{1},A_{1})w({\bf X}_{2},A_{2})\|{\bf X}_{1}-{\bf X}_{2}\|_{2}-\|{\bf
X}_{3}-{\bf X}_{4}\|_{2}\text{ and}$ $\displaystyle k_{A}(({\bf
X}_{1},A_{1}),\dots,({\bf X}_{4},A_{4}))={}$ $\displaystyle w({\bf
X}_{1},A_{1})|A_{1}-A_{3}|+w({\bf X}_{2},A_{2})|A_{2}-A_{4}|$
$\displaystyle-w({\bf X}_{1},A_{1})w({\bf
X}_{2},A_{2})|A_{1}-A_{2}|-|A_{3}-A_{4}|.$

The following assumptions are required for several Lemmas and Theorems
presented. These conditions amount to finite moment conditions on squares of
the Euclidean norms of ${\bf X}$, $A$ and their products with the weights
$w({\bf X},A)$.

  1. (A1)

$\mathbb{E}[k^{2}(({\bf X}_{1},A_{1}),\dots,({\bf X}_{6},A_{6}))]<\infty$

  2. (A2)

$\mathbb{E}[k_{\bf X}^{2}(({\bf X}_{1},A_{1}),\dots,({\bf
X}_{4},A_{4}))]<\infty$ and $\mathbb{E}[k_{A}^{2}(({\bf
X}_{1},A_{1}),\dots,({\bf X}_{4},A_{4}))]<\infty$

## References

  * (1)
  * Andersen et al. (2011) Andersen, M., Dahl, J., Liu, Z., Vandenberghe, L., Sra, S., Nowozin, S. and Wright, S. (2011), ‘Interior-point methods for large-scale cone programming’, Optimization for Machine Learning 5583.
  * Chattopadhyay et al. (2020) Chattopadhyay, A., Hase, C. H. and Zubizarreta, J. R. (2020), ‘Balancing vs modeling approaches to weighting in practice’, Statistics in Medicine 39(24), 3227–3254.
  * Colangelo and Lee (2020) Colangelo, K. and Lee, Y.-Y. (2020), ‘Double debiased machine learning nonparametric inference with continuous treatments’, arXiv preprint arXiv:2004.03036 .
  * Crump et al. (2009) Crump, R. K., Hotz, V. J., Imbens, G. W. and Mitnik, O. A. (2009), ‘Dealing with limited overlap in estimation of average treatment effects’, Biometrika 96(1), 187–199.
  * Dawid (1979) Dawid, A. P. (1979), ‘Some misleading arguments involving conditional independence’, Journal of the Royal Statistical Society: Series B (Methodological) 41(2), 249–252.
  * Díaz and van der Laan (2013) Díaz, I. and van der Laan, M. J. (2013), ‘Targeted data adaptive estimation of the causal dose–response curve’, Journal of Causal Inference 1(2), 171–192.
  * Fong et al. (2018) Fong, C., Hazlett, C. and Imai, K. (2018), ‘Covariate balancing propensity score for a continuous treatment: Application to the efficacy of political advertisements’, The Annals of Applied Statistics 12(1), 156–177.
  * Galvao and Wang (2015) Galvao, A. F. and Wang, L. (2015), ‘Uniformly semiparametric efficient estimation of treatment effects with a continuous treatment’, Journal of the American Statistical Association 110(512), 1528–1542.
  * Hill (2011) Hill, J. L. (2011), ‘Bayesian nonparametric modeling for causal inference’, Journal of Computational and Graphical Statistics 20(1), 217–240.
  * Hirano and Imbens (2004) Hirano, K. and Imbens, G. W. (2004), The Propensity Score with Continuous Treatments, John Wiley & Sons, Ltd, chapter 7, pp. 73–84.
  * Huling and Mak (2020) Huling, J. D. and Mak, S. (2020), ‘Energy balancing of covariate distributions’, arXiv preprint arXiv:2004.13962 .
  * Imai and Van Dyk (2004) Imai, K. and Van Dyk, D. A. (2004), ‘Causal inference with general treatment regimes: Generalizing the propensity score’, Journal of the American Statistical Association 99(467), 854–866.
  * Imbens (2004) Imbens, G. W. (2004), ‘Nonparametric estimation of average treatment effects under exogeneity: A review’, Review of Economics and Statistics 86(1), 4–29.
  * Johnson et al. (2016) Johnson, A. E., Pollard, T. J., Shen, L., Li-Wei, H. L., Feng, M., Ghassemi, M., Moody, B., Szolovits, P., Celi, L. A. and Mark, R. G. (2016), ‘Mimic-iii, a freely accessible critical care database’, Scientific data 3(1), 1–9.
  * Johnson et al. (2003) Johnson, E., Dominici, F., Griswold, M. and L. Zeger, S. (2003), ‘Disease cases and their medical costs attributable to smoking: an analysis of the national medical expenditure survey’, Journal of Econometrics 112(1), 135–151.
  * Kallus and Santacatterina (2019) Kallus, N. and Santacatterina, M. (2019), ‘Kernel optimal orthogonality weighting: A balancing approach to estimating effects of continuous treatments’, arXiv preprint arXiv:1910.11972 .
  * Kang and Schafer (2007) Kang, J. D. and Schafer, J. L. (2007), ‘Demystifying double robustness: A comparison of alternative strategies for estimating a population mean from incomplete data’, Statistical Science 22(4), 523–539.
  * Kennedy et al. (2017) Kennedy, E. H., Ma, Z., McHugh, M. D. and Small, D. S. (2017), ‘Non-parametric methods for doubly robust estimation of continuous treatment effects’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 79(4), 1229–1245.
  * Kish (1965) Kish, L. (1965), Survey sampling, Vol. 26, John Wiley & Sons.
  * Mak and Joseph (2018) Mak, S. and Joseph, V. R. (2018), ‘Support points’, The Annals of Statistics 46(6A), 2562–2592.
  * Martinet (2020) Martinet, G. (2020), ‘A balancing weight framework for estimating the causal effect of general treatments’, arXiv preprint arXiv:2002.11276 .
  * Naimi et al. (2014) Naimi, A. I., Moodie, E. E., Auger, N. and Kaufman, J. S. (2014), ‘Constructing inverse probability weights for continuous exposures: a comparison of methods’, Epidemiology pp. 292–299.
  * Neto et al. (2018) Neto, A. S., Deliberato, R. O., Johnson, A. E., Bos, L. D., Amorim, P., Pereira, S. M., Cazati, D. C., Cordioli, R. L., Correa, T. D., Pollard, T. J. et al. (2018), ‘Mechanical power of ventilation is associated with mortality in critically ill patients: an analysis of patients in two observational cohorts’, Intensive Care Medicine 44(11), 1914–1922.
  * Papazian et al. (2019) Papazian, L., Aubron, C., Brochard, L., Chiche, J.-D., Combes, A., Dreyfuss, D., Forel, J.-M., Guérin, C., Jaber, S., Mekontso-Dessap, A. et al. (2019), ‘Formal guidelines: management of acute respiratory distress syndrome’, Annals of Intensive Care 9(1), 1–18.
  * Robins et al. (2000) Robins, J. M., Hernan, M. A. and Brumback, B. (2000), ‘Marginal structural models and causal inference in epidemiology’, Epidemiology 11, 550–560.
  * Rosenbaum and Rubin (1983) Rosenbaum, P. R. and Rubin, D. B. (1983), ‘The central role of the propensity score in observational studies for causal effects’, Biometrika 70(1), 41–55.
  * Siegel et al. (2019) Siegel, M. D., Hyzy, R. C. et al. (2019), ‘Ventilator management strategies for adults with acute respiratory distress syndrome’.
  * Stellato et al. (2020) Stellato, B., Banjac, G., Goulart, P., Bemporad, A. and Boyd, S. (2020), ‘Osqp: An operator splitting solver for quadratic programs’, Mathematical Programming Computation pp. 1–36.
  * Székely and Rizzo (2013) Székely, G. J. and Rizzo, M. L. (2013), ‘Energy statistics: A class of statistics based on distances’, Journal of Statistical Planning and Inference 143(8), 1249–1272.
  * Székely et al. (2007) Székely, G. J., Rizzo, M. L. and Bakirov, N. K. (2007), ‘Measuring and testing dependence by correlation of distances’, The Annals of Statistics 35(6), 2769–2794.
  * Tübbicke (2020) Tübbicke, S. (2020), ‘Entropy balancing for continuous treatments’, arXiv preprint arXiv:2001.06281 .
  * van der Laan and Robins (2003) van der Laan, M. J. and Robins, J. M. (2003), Unified methods for censored longitudinal data and causality, Springer Science & Business Media.
  * Vegetabile et al. (2021) Vegetabile, B. G., Griffin, B. A., Coffman, D. L., Robbins, M. W., Cefalu, M. and McCaffrey, D. F. (2021), ‘Nonparametric estimation of population average dose-response curves using entropy balancing weights for continuous exposures’, Health Services and Outcomes Research Methodology 21(1), 69–110.
  * Wang and Wahba (1995) Wang, Y. and Wahba, G. (1995), ‘Bootstrap confidence intervals for smoothing splines and their comparison to bayesian confidence intervals’, Journal of Statistical Computation and Simulation 51(2-4), 263–279.
  * Yiu and Su (2018) Yiu, S. and Su, L. (2018), ‘Covariate association eliminating weights: a unified weighting framework for causal effect estimation’, Biometrika 105(3), 709–722.
  * Zhu et al. (2015) Zhu, Y., Coffman, D. L. and Ghosh, D. (2015), ‘A boosting algorithm for estimating generalized propensity scores with continuous treatments’, Journal of Causal Inference 3(1), 25–40.
  * Zubizarreta (2015) Zubizarreta, J. R. (2015), ‘Stable weights that balance covariates for estimation with incomplete outcome data’, Journal of the American Statistical Association 110(511), 910–922.