# Causal Bias Quantification for Continuous Treatments

Gianluca Detommaso  Michael Brückner  Philip Schulz  Victor Chernozhukov

###### Abstract

We extend the definition of the marginal causal effect to the continuous
treatment setting and develop a novel characterization of causal bias in the
framework of structural causal models. We prove that our derived bias
expression is zero if, and only if, the causal effect is identifiable via
covariate adjustment. We show that under some restrictions on the structural
equations, the causal bias can be estimated efficiently and allows for causal
regularization of predictive probabilistic models. We demonstrate the
effectiveness of our method for causal bias quantification in various settings
where (not) controlling for certain covariates would introduce causal bias.

Causal Inference

##  1 Introduction

In probabilistic modeling we are often concerned with estimating the posterior
predictive distribution $p(Y|X=x,\mathcal{D})$ from observed data
$\mathcal{D}$, generated by some unknown data generating mechanism (Bishop,
2006; Koller & Friedman, 2009). The posterior predictive captures our beliefs
about the association of some input variable $X$ and the corresponding outcome
variable $Y$ when both were generated by the same mechanism as the training
data. However, in many practical applications, new inputs might not be
generated by the exact same mechanism; we may face selection bias; or we aim
for a fair model explaining outcomes primarily based on cause. In these cases,
we may consider a causal model and estimate the counterfactual distribution
(posterior interventional distribution), $p(Y|X\leftarrow x,\mathcal{D})$
(Lauritzen, 2001; Pearl, 2009a, b), which states the posterior probability of
an outcome if we intervene in the data generating mechanism and set treatment
variable (input) $X$ to $x$, also written $do(X=x)$.

To illustrate the difference, consider watching people (not) wearing a
raincoat and we are interested in predicting street conditions. In a purely
predictive setting, we would learn from past data that observing raincoats is
correlated with wet streets and predict street conditions accordingly.
However, such a model may fail when applied to a different population of
people. For instance, in windy regions, people may wear raincoats even on
sunny days; in warm regions, people may prefer umbrellas. In contrast, the
posterior interventional distribution would tell us that putting on a coat
does not affect street conditions and our prediction would purely rely on the
marginal probability of observing wet streets.

In many applications, the situation is more nuanced and we are facing the
problem whether to include certain covariates in the model to improve
predictive performance or whether to exclude them to not incur into bias. The
problem is well-understood if we aim to avoid any causal bias; however, what
if we are willing to accept some bias in order to improve predictive
performance?

Starting from a given causal probabilistic model defined in terms of
structural equations, we are interested in quantifying the causal bias of a
predictive model which, together with its posterior probability, allows us to
balance predictive performance with causal alignment. We focus on the setting
of continuous treatment variables, which will enable us to state an analytic
expression of the causal bias.

###  1.1 Related Work

Large parts of literature on causal inference focus on causal identifiability,
that is, identifying conditions under which the causal relationship can be
unbiasedly estimated from the observed data (Spirtes et al., 2000; Shpitser &
Pearl, 2006; Pearl, 2009c), and causal estimation of models when these
conditions are satisfied.

If we are not willing to make any additional (parametric) modeling
assumptions, causal identifiability requires treatments from a finite domain
(Angrist & Pischke, 2008; Hernán & Robins, 2020). In particular, for binary
treatments, we might only be interested in the causal effect expressed in
terms of the average treatment effect of the treated (ATET) or the average
treatment effect (ATE) (Imbens & Rubin, 2015), which are the expected
difference between the (potential) outcomes under the two alternative
treatments for the treated and the whole population, respectively.

In situations where different posterior interventional distributions are
compatible with the observed data and the assumed causal relations, the causal
effect is non-identifiable without making additional modeling assumptions, for
instance, on the functional relationship of treatments and outcomes; this is
in particular true for continuous treatments (Rosenbaum et al., 2010; Silva,
2016; Hines et al., 2021). Modeling the treatment assignment mechanism allows
for partial identification of the causal effect through inverse propensity
weighting (Imbens, 2000; Hirano & Imbens, 2004; Imai & Van Dyk, 2004; Wu et
al., 2018; Austin, 2019), modeling the outcome mechanism allows for regression
methods (Chipman et al., 2010; Hill, 2011), and doubly robust estimation can
be used to combine both (Kennedy et al., 2017; Chernozhukov et al., 2018).
More generally, structural causal models (Pearl, 2009c) allow us to describe
the full data generating mechanism.

Similarly, estimating the average treatment effect for continuous treatments
either requires some discretization or assuming a certain functional form
called dose-response function in case of bounded univariate treatments
(Altshuler, 1981; Gill & Robins, 2001; Wang, 2015; Galagate, 2016).

However, most existing works on continuous treatments assume that the causal
effect is identifiable and, hence, the causal bias is zero. In contrast, we
consider the scenario where the causal effect might be non-identifiable and
develop a characterization of the bias that can be used to evaluate how far a
model is from being causal. While different in methodology, a similar purpose
can be found in (Tran et al., 2016; Gain & Shpitser, 2018). In addition, it
can be used to perform sensitivity analysis in the presence of unobserved
confounders in a completely non-linear framework, similarly to what is
proposed in (Cinelli & Hazlett, 2020) for the linear case. Our proposed
criterion for identifiability via covariate adjustment is also related to the
adjustment criterion in (Shpitser et al., 2010), which was shown by the
authors to imply and be implied by conditional ignorability (Rosenbaum &
Rubin, 1983; Little & Rubin, 2019). Our criterion is useful to build the link
between identifiability and causal bias being zero.

###  1.2 Preliminaries

We assume a probabilistic model of the data generating mechanism, that is, a
joint distribution over some random variables
$\mathcal{V}=\\{V_{1},\ldots,V_{n}\\}$ and a causal directed acyclic graph
(DAG) $\mathcal{G}$ encoding the direct causal relationship between them
(Lauritzen, 2001). We denote parents, ancestors, and descendants in this graph
by $\text{Pa}(V_{i})\subset\mathcal{V}$, $\text{An}(V_{i})\subset\mathcal{V}$,
and $\text{De}(V_{i})\subset\mathcal{V}$, respectively.

We further assume that the probabilistic model and causal graph are stated in
terms of a canonical structural causal model (SCM) (Pearl, 2009c) consisting
of a set of independent random variables
$\mathcal{U}=\\{U_{V_{1}},\ldots,U_{V_{n}}\\}$ with some probability densities
$p(U_{V_{i}})$, and a set of functions
$\mathcal{F}=\\{f_{V_{1}},\ldots,f_{V_{n}}\\}$ such that
$V_{i}:=f_{V_{i}}(\text{Pa}(V_{i}),U_{V_{i}})\ \forall i$. The exogenous
variables $U_{V_{i}}$ are latent and cannot be observed, whereas the
endogenous variables $V_{i}$ may or may not be observed.

Technically, an SCM defines a joint probability distribution over
$\mathcal{U}$ and $\mathcal{V}$ satisfying the global Markov property, where
marginalizing over $\mathcal{U}$ yields a joint observational distribution
over all endogenous variables $\mathcal{V}$. We denote the graph corresponding
to an SCM by $\mathcal{G}^{+}$ and denote parents, ancestors, and descendants
in $\mathcal{G}^{+}$ accordingly by
$\text{Pa}^{+}(\cdot)\subset\mathcal{V}\cup\mathcal{U}$,
$\text{An}^{+}(\cdot)\subset\mathcal{V}\cup\mathcal{U}$, and
$\text{De}^{+}(\cdot)\subset\mathcal{V}\cup\mathcal{U}$. Figure 1(a) shows an
example of a causal graph $\mathcal{G}$ and Figure 1(b) shows a graph
$\mathcal{G}^{+}$ of an SCM that aligns with $\mathcal{G}$.

$V_{1}$$V_{2}$$V_{3}$$V_{4}$$V_{5}$ (a) Example of a causal graph
$\mathcal{G}$.

$V_{1}$$V_{2}$$V_{3}$$V_{4}$$V_{5}$$U_{V_{1}}$$U_{V_{2}}$$U_{V_{3}}$$U_{V_{4}}$$U_{V_{5}}$
(b) Example of reparameterized causal graph $\mathcal{G}^{+}$.

Figure 1: Example of a causal graph (a) and a (reparameterized) causal graph
corresponding to an SCM that agrees with $\mathcal{G}$. Grey and white circles
denote observed and unobserved variables, respectively. Stochastic
dependencies are encoded through solid arrows whereas dashed arrows highlight
that the transitions are deterministic given all parents.

Among the endogenous variables, we assume a single, potentially multivariate,
treatment variable $X\in\mathcal{V}$ with $\mathcal{X}:=\\{X\\}$, and a single
outcome variable $Y\in\mathcal{V}$ with $\mathcal{Y}:=\\{Y\\}$. We further
denote the remaining set of observed endogenous variables by
$\mathcal{O}\subseteq\mathcal{V}\setminus(\mathcal{X}\cup\mathcal{Y})$, while
all other endogenous variables
$\mathcal{L}:=\mathcal{V}\setminus(\mathcal{X}\cup\mathcal{Y}\cup\mathcal{O})$
are latent. Among the observed variables $\mathcal{O}$, we denote non-
descendants of $X$ by
$\underline{\mathcal{O}}:=\mathcal{O}\setminus\text{De}(X)$, similarly to the
notation in mutilated graphs (Pearl, 2009b). We denote the corresponding set
of exogenous variables by $\mathcal{U}_{\mathcal{O}}$,
$\mathcal{U}_{\underline{\mathcal{O}}}$, and $\mathcal{U}_{\mathcal{L}}$,
respectively.

We use lowercase characters to denote realizations of random variables, for
instance, $v_{i}$ for a realization of $V_{i}$, and we may use
$p(\cdot|v_{i})$ as a short form for the conditional probability density (or
mass) $p(\cdot|V_{i}=v_{i})$. We use capital letters without indices to denote
concatenated vectors of random variables, for instance, $O:=[O_{j}]_{j=1}^{m}$
with $\mathcal{O}:=\\{O_{1},\ldots,O_{m}\\}$, and
$U_{\mathcal{O}}:=[U_{O_{j}}]_{j=1}^{m}$. Likewise, realizations without
indices refer to vectors of realizations like $o:=[o_{j}]_{j=1}^{m}$, and
$f_{O}$ denotes the concatenation of all $f_{O_{j}}$, which is a vector-valued
function with combined set of arguments.

For a generic function $h$ and argument(s) $x$, the Jacobian of $h$ with
respect to $x$ is denoted by $\nabla_{x}h$ if it exists. Note that if $h$ is a
scalar function, $\nabla_{x}h$ denotes a transposed gradient. We use
$\mathbb{E}$, Cov and Var for expectation, covariance, and variance. We may
use subscripts to specify the underlying probability density, e.g.
$\mathbb{E}_{Y|x,o}$ means that the expectation is taken with respect to
$p(Y|x,o)$.

##  2 Identification, quantification and estimation

In this section we (i) introduce a novel criterion for identifiability of
causal effects via covariate adjusts; (ii) we derive closed-form expressions
for causal effect and causal bias that can be computed in closed-form; and
(iii) show how the latter can be estimated via standard probabilistic
inference methods.

###  2.1 A characterization of identifiability via covariate adjustment

This section offers a complete criterion for identifiability of causal effects
(Pearl, 2009b) via covariate adjustment (Shpitser et al., 2010) in the
framework of structural causal models. The underlying idea is to first
construct a minimal set of exogenous variables,
$\mathcal{U}^{C}\subseteq\mathcal{U}$, such that, together with variables in
$\mathcal{X}$ and $\underline{\mathcal{O}}$, this set fully characterizes the
randomness of outcome variable $Y$. To this end, let us again consider the
causal graph in Figure 1(a) as a running example.

######  Example 1.

We assume the causal graph in Figure 1(a) and $X=V_{2}$, $Y=V_{4}$,
$\mathcal{O}=\\{V_{3},V_{5}\\}$, and $\mathcal{L}=\\{V_{1}\\}$. In this case,
it is well understood that the unobserved variable $V_{1}$ and the observed
variables $V_{3}$ and $V_{5}$ each compromise identifiability of the causal
effect from $X$ to $Y$. Latent variable $V_{1}$ is a hidden confounder
introducing confounding bias, observing mediator $V_{3}$ introduces
overcontrol bias, and observing collider $V_{5}$ introduces (endogenous)
selection bias (Elwert, 2013). Figure 1(b) is the causal graph
$\mathcal{G}^{+}$ of an SCM aligning with $\mathcal{G}$. We notice that

$U_{V_{1}}\not\perp\\!\\!\\!\perp X,\quad U_{V_{3}}\not\perp\\!\\!\\!\perp
X|V_{3},\quad U_{Y}\not\perp\\!\\!\\!\perp X|V_{5},$

which we will see in Theorem 2.4 to be the very reasons why the causal effect
is not identifiable through covariate adjustment.

######  Definition 2.1.

For any set of treatment variables $\mathcal{X}$ and outcome variable $Y$, the
set of exogenous random variables which are ancestors of $Y$ in
$\mathcal{G}^{+}$ and that are not independent of $Y$ given $\mathcal{X}$ and
$\underline{\mathcal{O}}$ is called the causal exogenous set:

$\mathcal{U}^{C}:=\\{U_{V_{i}}\in\mathcal{U}\cap\text{An}^{+}(Y):U_{V_{i}}\not\perp\\!\\!\\!\perp
Y|\mathcal{X},\underline{\mathcal{O}}\\}.$ (1)

Intuitively, this set contains all exogenous variables which we would need to
control for to translate the observational setting into an experimental
setting. In Example 1 this set includes random variables $U_{V_{1}}$,
$U_{V_{3}}$, and $U_{Y}$. While we can never observe exogenous variables, we
might be able to obtain a posterior distribution over them by making (weak)
prior distributional assumptions. If this posterior distribution is
insensitive to the treatment variable, the causal effect can be (partially)
identified. In the following, we state this observation more formally.

######  Assumption 2.2.

We assume that the given SCM is defined through continuous exogenous variables
$U_{V_{i}}$ where the probability density exists and is strictly positive for
all possible values of $U_{V_{i}}$.

######  Assumption 2.3.

We assume that the given SCM is faithful to the causal graph
$\mathcal{G}^{+}$, that is, any independences that hold in the model must hold
in the causal graph.

######  Theorem 2.4.

Given an SCM that satisfies Assumption 2.2 with causal exogenous set
$\mathcal{U}^{C}$ for treatment variable $X$, outcome variable $Y$, and a set
of observations $\mathcal{O}$, the causal effect of $X$ on $Y$ can be
identified by covariate adjustment if, and only if,

$U_{V_{i}}\perp\\!\\!\\!\perp X|\mathcal{O}\ \ \ \forall\
U_{V_{i}}\in\mathcal{U}^{C}$ (2)

###### Proof.

See Appendix A.1. ∎

######  Corollary 2.5.

Given an SCM that satisfies Assumptions 2.2-2.3 where the causal effect of $X$
on $Y$ can be identified by covariate adjustment and which is differentiable
in $X$, then

$\left.\nabla_{X}\log p(U^{C}|X,O)\right|_{X=x}=0\ \ \ \forall\ x,$ (3)

where $U^{C}$ denotes the vector of all endogenous variables in
$\mathcal{U}^{C}$ and $O$ denotes the vector of observations.

###### Proof.

According to Theorem 2.4, treatment variable $X$ is independent of all
elements of the causal exogenous set given the observed variables,
$p(U^{C}|X,O)=p(U^{C}|O)$. If the density defined through the SCM exists and
is differentiable w.r.t. $X$, the partial gradient $\nabla_{X}\log
p(U^{C}|X,O)$ is zero everywhere. ∎

###  2.2 Quantification of average partial effect and bias

Corollary 2.5 gives rise to a measure of causal bias which is zero if the
causal bias is zero. To formalize this, let us introduce the notation of
partial evaluation, which is related to the potential outcome notation and
interventions in structural causal models.

######  Definition 2.6.

Given a function $f$ with arguments
$\mathcal{A}=\\{V_{i_{1}},\ldots,V_{i_{m}}\\}\subseteq\mathcal{V}$ and some
set $\mathcal{B}=\\{V_{j_{1}},\ldots,V_{j_{k}}\\}\subseteq\mathcal{V}$, we
define $f^{b}$ to be the partial evaluation of $f$ where all arguments
$\mathcal{A}\cap\mathcal{B}$ have been assigned to the corresponding value
$b=[b_{1},\ldots,b_{k}]$, that is

$f^{\mathcal{B}\leftarrow
b}(\mathcal{A}\setminus\mathcal{B}):=f(V_{i_{1}},\dots,V_{i_{m}})\big{|}_{V_{j_{1}}=b_{1},\ldots,V_{j_{k}}=b_{k}}.$
(4)

We may use the short form $f^{b}$ for $f^{\mathcal{B}\leftarrow
b}(\mathcal{A}\setminus\mathcal{B})$ if the corresponding set of assigned
random variables $\mathcal{B}$ is clear from the context, for example, all
observed variables $\mathcal{O}$. Likewise, if sets $\mathcal{B}$ and
$\mathcal{C}$ are disjoint sets of random variables, we use $f^{b,c}$ to
denote simultaneously assignments $\mathcal{B}\leftarrow b$ and
$\mathcal{C}\leftarrow c$.

######  Example 2.

Consider again the causal graph in Figure 1(b) with $X=V_{2}$ with realization
$x$, $\mathcal{O}=\\{V_{3},V_{5}\\}$ with realization $o=[o_{1},o_{2}]$,
$Y=V_{4}$, and $\mathcal{L}=\\{V_{1}\\}$. Here, for instance, we have

$\displaystyle f_{Y}^{x}$ $\displaystyle=$ $\displaystyle
f_{Y}(V_{1},x,f_{V_{3}}(x,U_{V_{3}}),U_{Y}),$ $\displaystyle f_{V_{5}}^{x,o}$
$\displaystyle=$ $\displaystyle
f_{V_{5}}(x,f_{Y}(V_{1},x,o_{1},U_{Y}),U_{V_{5}}).$

If derivatives are well-defined, this for example implies
$\nabla_{U_{V_{3}}}f_{V_{5}}^{x,o}=0$, while in general
$\nabla_{U_{V_{3}}}f_{V_{3}}\neq 0$.

######  Proposition 2.7.

Given an SCM that satisfies Assumptions 2.2-2.3 with a vector of causal
exogenous variables $U^{C}$ for treatment variable $X$, outcome variable $Y$,
and observations $O$. Suppose that $f_{Y}$ and $p$ are differentiable with
respect to $X$ and satisfy sufficient regularity assumptions. Then

$\displaystyle\underbrace{\nabla_{x}\mathbb{E}_{Y|x,o}[Y]}_{\mathcal{A}(x,o)}$
$\displaystyle\\!\\!=\\!\\!$
$\displaystyle\underbrace{\mathbb{E}_{U^{C}|x,o}\left[\nabla_{x}f_{Y}^{x,\underline{o}}\right]}_{\mathcal{C}(x,o)}\
+$ (5)
$\displaystyle\underbrace{\mathbb{E}_{U^{C}|x,o}\left[f_{Y}^{x,\underline{o}}\nabla_{x}\\!\log
p(U^{C}|x,o)\right]}_{\mathcal{B}(x,o)}\\!.\ \ \ \ $

If the causal effect of $X$ on $Y$ can be identified by covariate adjustment
then $\mathcal{B}(x,o)=0$.

###### Proof.

See Appendix A.2. ∎

In the above proposition, expression $\mathcal{A}(x,o)$ states the partial
derivative of the average observational outcome, $\mathcal{C}(x,o)$ states
average partial effect on the treated (APET) for an infinitesimal change in
the treatment, and $\mathcal{B}(x,o)$ is the causal bias.

######  Remark 2.8.

The decomposition in (5) has a well-known discrete counterpart for some
reference treatment $x$ (control) and proposal $x^{\prime}:=x+h$ (treated),

$\displaystyle\underbrace{\mathbb{E}[Y|x^{\prime}]-\mathbb{E}[Y|x]}_{\text{association}}$
$\displaystyle\\!\\!=\\!\\!$
$\displaystyle\underbrace{\mathbb{E}[Y^{x^{\prime}}-Y^{x}|x^{\prime}]}_{\text{ATET}}\
+$
$\displaystyle\underbrace{\mathbb{E}[Y^{x}|x^{\prime}]-\mathbb{E}[Y^{x}|x]}_{\text{bias}}\,,\
$

where association is defined as the contrast in average outcomes by observed
treatment value, and ATET refers to the average treatment effect on the
treated. If we condition on $O=o$, divide both sides by $h$ and take the limit
for $h\to 0$, we recover (5).

The expressions for $\mathcal{C}(x,o)$ and $\mathcal{B}(x,o)$ in (5) highlight
the connection with Theorem 2.4. However, in practice we do not need to find
the causal exogenous set $\mathcal{U}^{C}$ in order to compute them. Indeed,
because $\nabla_{x}f_{Y}^{x,\underline{o}}$ does not depend on
$\mathcal{U}\setminus\mathcal{U}^{C}$ by definition of $\mathcal{U}^{C}$, via
marginalization we can write

$\mathcal{C}(x,o)=\mathbb{E}_{U|x,o}\left[\nabla_{x}f_{Y}^{x}\right],$ (6)

where $U$ denotes the vector of all exogenous variables and we are free to
remove the assignment of $\underline{o}$. Regarding $\mathcal{B}$, the main
challenge of the expression in (5) is to compute the term $\nabla_{x}\log
p(U^{C}|x,o)$, since $x$ appears in the normalization constant. Theorem 2.9
provides an elegant formulation of the bias that is easier to compute in
practice. Under some additional assumptions Theorem 2.10 provides an
alternative formulation of $\mathcal{B}$ that is automatically differentiable
and does not depend on $U^{C}$.

######  Theorem 2.9.

Under the assumptions of Proposition 2.7, the bias in (5) can be rewritten as

$\mathcal{B}(x,o)=\text{Cov}_{U^{C}|x,o}\left(f_{Y}^{x,\underline{o}},\nabla_{x}\log
p(U^{C},x,o)\right),$ (7)

where $p(U^{C},x,o)$ denotes the joint density over $U^{C}$, $X$, and $O$
(partially) evaluated for $x$ and $o$.

###### Proof.

See Appendix A.3. ∎

The expression of the bias in Theorem 2.9 is a clear improvement over (5) as
the joint $p(U^{C},X,O)$ does not include a normalization constant involving
$X$. If no mediators between $X$ and $Y$ are observed (e.g. $V_{3}$ in Example
1), (7) can already provide a practical way to estimate the bias. Otherwise it
may hide a difficulty: If the factorization of the joint includes a Dirac
delta distribution depending on $X$, computing its gradient with respect to
$X$ requires some extra analysis.

For instance, consider Example 1 where the joint contains conditional
$p(V_{3}|X,U_{V_{3}})$ as $V_{3}\in\mathcal{O}$. Since this conditional is a
Dirac delta distribution, how can we then compute the gradient of $\log
p(U^{C},X,O)$ w.r.t. $X$? The following theorem further manipulates the bias
to remove this obstacle.

######  Theorem 2.10.

Suppose that $f_{Y}$ is differentiable in $X$ and all $V_{i}\in\mathcal{O}$,
and that corresponding functions $f_{X}$ and $f_{V_{i}}$ are differentiable
and invertible in $U_{X}$ and $U_{V_{i}}$, respectively. Then the bias in (5)
can be rewritten as

$\displaystyle\mathcal{B}(x,o)\
=\\!\\!\sum_{V_{i}\in\mathcal{X}\cup\mathcal{O}}\\!\\!\\!\\!\\!\\!\mathbb{E}_{U|x,o}\Big{[}$
$\displaystyle\Big{(}\nabla_{U_{V_{i}}}f_{Y}^{x}\\!+\\!\left(f_{Y}^{x}-y^{x}\right)\nabla_{U_{V_{i}}}\\!\log
p(U_{V_{i}})\Big{)}$ $\displaystyle\
(\nabla_{U_{V_{i}}}f_{V_{i}}^{x,o})^{-1}\,\nabla_{x}(f_{V_{i}}^{x,o}-v_{i})\Big{]}$
(8)

with potential outcome $y^{x}:=\mathbb{E}_{U|x,o}[f_{Y}^{x}]$ for treatment
$x$.

###### Proof.

See Appendix A.4. ∎

######  Remark 2.11.

Note that (i) all terms in the expectation of (2.10) can be efficiently
computed via automatic differentiation. (ii) Bias $\mathcal{B}(x,o)$ can be
decomposed in the contributions given by $X$ (potential confounding bias) and
by each variable $V_{i}\in\mathcal{O}$ (potential overcontrol or selection
bias). (iii) If $U_{V_{i}}$ is modeled as a (multivariate) standard Gaussian,
then $\nabla_{U_{V_{i}}}\log p(U_{V_{i}})=-U_{V_{i}}$. Finally, (iv) the
inverse Jacobian $(\nabla_{U_{V_{i}}}f_{V_{i}}^{x,o})^{-1}$ can be efficiently
computed if $V_{i}$ is low-dimensional or if $f_{V_{i}}$ has, for example, a
normalizing flow structure (Papamakarios et al., 2019).

###  2.3 Estimation of average partial effect and bias

Sections 2.1 and 2.2 provide closed-form expressions for average partial
effect on the treated and the causal bias, namely (6) and (2.10). Their
numerical estimation reduces to a mere statistical inference task: the
approximation of an expectation with respect to the probability density
$p(U|x,o)$. We further notice that

$p(U|x,o)=p(U_{\mathcal{L}}|x,o)\prod_{V_{i}\in\mathcal{X}\cup\mathcal{O}}\delta(U_{V_{i}}-u_{V_{i}}^{*}),$
(9)

where $u_{V_{i}}^{*}$ are uniquely determined by the assumptions in Theorem
2.10 that $f_{V_{i}}$ are invertible in $U_{V_{i}}$ for all
$V_{i}\in\mathcal{X}\cup\mathcal{O}$ given $U_{\mathcal{L}}$. Hence, we can
restrict our attention to infer $p(U_{\mathcal{L}}|x,o)$. In some cases, such
a posterior might be available in closed-form. Otherwise we may resort to
sampling methods such as MCMC (Brooks et al., 2011), SMC (Doucet et al.,
2001), and importance sampling (Neal, 2001), or variational methods like
Laplace approximation (Shun & McCullagh, 1995), variational Bayes (Jordan et
al., 1999; Blei et al., 2017), and normalizing flows (Rezende & Mohamed, 2015;
Papamakarios et al., 2019). Let us quickly review Laplace approximation and
importance sampling in our framework, as we will use them in Sections 3.2 and
3.3, respectively.

Laplace approximation. The Laplace approximation is a Gaussian approximation
of $p(U_{\mathcal{L}}|x,o)$ around its mode. In order to compute the mode, we
first notice that via automatic differentiation we have immediate access to
$\nabla_{U_{V_{i}}}\log p(U_{\mathcal{L}},x,o)$ for all
$U_{V_{i}}\in\mathcal{U}_{\mathcal{L}}$, which we can employ in any gradient-
based optimization to compute a maximum-a-posteriori (MAP) estimate (Murphy,
2012). Let us denote the latter by $u_{\mathcal{L}}^{\scriptsize\text{MAP}}$.
Furthermore, automatic differentiation also gives access to the Hessian with
elements (sub-matrices) $\nabla_{U_{V_{i}}}\nabla_{U_{V_{j}}}\log
p(U_{\mathcal{L}},x,o)$ or to any scalable approximation of it (e.g. diagonal
approximation). Denote this Hessian evaluated at the MAP estimate by
$H(u_{\mathcal{L}}^{\scriptsize\text{MAP}})$. Whenever
$p(U_{\mathcal{L}}|x,o)$ has a well-defined MAP estimate,
$-H(u_{\mathcal{L}}^{\scriptsize\text{MAP}})$ is positive-definite and we can
approximate the posterior by a (multivariate) Gaussian,

$p(U_{\mathcal{L}}|x,o)\approx\mathcal{N}\\!\left(u_{\mathcal{L}}^{\scriptsize\text{MAP}},-H^{-1}(u_{\mathcal{L}}^{\scriptsize\text{MAP}})\right)\\!(U_{\mathcal{L}}).$
(10)

We can now generate a sample $u_{\mathcal{L}}^{(i)}$ from the approximated
posterior, compute the remaining exogenous variables via the respective
inverse maps of $f_{V_{i}}^{x,o}$, and produce a Monte Carlo estimate of
$\mathcal{C}(x,o)$ and $\mathcal{B}(x,o)$ using (6) and (2.10), respectively.

Importance sampling. Importance sampling is a way to modify the density
underlying an expectation. In causal inference, it is used for inverse
propensity weighting (IPW) (Imbens, 2000; Hirano & Imbens, 2004; Imai & Van
Dyk, 2004; Wu et al., 2018; Austin, 2019). Consider any treatable density
$q(U_{\mathcal{L}})>0$ of a sampling distribution, and some function
$h(u_{\mathcal{L}},u_{X},u_{\mathcal{O}})$, which is integrable in all
arguments. We have

$\displaystyle\mathbb{E}_{U|x,o}[h(U_{\mathcal{L}},\\!\\!\\!\\!\\!\\!\\!\\!$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!U_{X},U_{\mathcal{O}})]$ (11)
$\displaystyle=$
$\displaystyle\mathbb{E}_{U_{\mathcal{L}}|x,o}\left[h(U_{\mathcal{L}},u_{X}^{*},u_{\mathcal{O}}^{*})\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{U_{\mathcal{L}}}\left[h(U_{\mathcal{L}},u_{X}^{*},u_{\mathcal{O}}^{*})w(U_{\mathcal{L}}|x,o)\right],\
\ \ \ \ \ $

where $\mathbb{E}_{U_{\mathcal{L}}}$ denotes the expectation with respect to
density $q(\mathcal{U}_{\mathcal{L}})$ and normalized importance weights

$w(U_{\mathcal{L}}|x,o)\propto\frac{p(U_{\mathcal{L}},x,o)}{q(U_{\mathcal{L}})},$
(12)

with $\mathbb{E}_{U_{\mathcal{L}}}\left[w(U_{\mathcal{L}}|x,o)\right]=1$. Then
we can first sample $u_{\mathcal{L}}^{(i)}$ from $q(U_{\mathcal{L}})$ for
$i=1,\ldots,N$, compute the unnormalized weights
$\tilde{w}^{(i)}=p(u_{\mathcal{L}}^{(i)},x,o)\,/\,q(u_{\mathcal{L}}^{(i)})$,
and normalize the importance weights to sum up to $N$. Finally, using these
same samples, we can compose a Monte Carlo estimator using (11), and estimate
$\mathcal{C}(x,o)$ and $\mathcal{B}(x,o)$ by applying this technique for $h$
being the respective expressions in (6) and (2.10). Notice that if we take
$q(U_{\mathcal{L}})$ to be the prior $p(U_{\mathcal{L}})$, the unnormalized
weights simplify to be the likelihood $p(x,o|u_{\mathcal{L}}^{(i)})$.

##  3 Experiments

$V_{1}$$X$$Y$$U_{V_{1}}$$U_{X}$$U_{Y}$

$\displaystyle V_{1}$ $\displaystyle=U_{V_{1}}$ $\displaystyle X$
$\displaystyle=\alpha V_{1}+U_{X}$ $\displaystyle Y$ $\displaystyle=\beta
X+\gamma V_{1}+U_{Y}$

Figure 2: Confounding bias

$X$$V_{1}$$Y$$U_{X}$$U_{V_{1}}$$U_{Y}$

$\displaystyle X$ $\displaystyle=U_{X}$ $\displaystyle V_{1}$
$\displaystyle=\alpha X+U_{V_{1}}$ $\displaystyle Y$ $\displaystyle=\beta
X+\gamma V_{1}+U_{Y}$

Figure 3: Overcontrol bias

$X$$Y$$V_{1}$$U_{V_{1}}$$U_{X}$$U_{Y}$

$\displaystyle X$ $\displaystyle=U_{X}$ $\displaystyle Y$
$\displaystyle=\alpha X+U_{Y}$ $\displaystyle V_{1}$ $\displaystyle=\beta
X+\gamma Y+U_{V_{1}}$

Figure 4: End. selection bias

In this section we show several experiments that validate the theory developed
in Section 2, and illustrate interesting applications. In Section 3.1, we
provide closed-form expressions for all possible graphical structure inducing
causal bias, and linear structure equation models. in Section 3.2, we show how
quantifying the bias can be helpful in situation where full identifiability
cannot be achieved due to missing data. Finally, in Section 3.3 we study the
application of our framework to a simulated medical study.

###  3.1 Confounding, overcontrol and selection bias

First, we study simple linear models exhibiting the three possible types of
bias: confounding, overcontrol and (endogenous) selection (Elwert, 2013).
Studying linear models is useful because we can compute $\mathcal{A}$,
$\mathcal{C}$ and $\mathcal{B}$ in closed-form, checking our results and
building up intuition. For all models, we will consider $\alpha,\beta,\gamma$
and $\delta$ to be scalar parameters and
$U_{V_{1}},U_{X},U_{Y}\sim\mathcal{N}(0,1)$. For detailed proofs of the
following results, see Appendix B.

Confounding bias. Given the model in Figure 4, we can compute

$\mathcal{A}(x,o)\,=\,\beta+\frac{\gamma\alpha}{1+\alpha^{2}},\quad\mathcal{C}(x,o)\,=\,\beta,\quad\mathcal{B}(x,o)\,=\,\frac{\gamma\alpha}{1+\alpha^{2}}.$

Notice that because the model is linear, $\mathcal{C}(x,o)$ is independent of
$p(U|x,o)$. If $\gamma\neq 0$ then $U_{V_{1}}\in\mathcal{U}^{C}$, and if
$\alpha\neq 0$ then $U_{V_{1}}\not\perp\\!\\!\\!\perp X$. Thus (2) is violated
and $\mathcal{B}(x,o)\neq 0$.

Overcontrol bias. Given the model in Figure 4, we can compute

$\mathcal{A}(x,o)\,=\,\beta,\quad\mathcal{C}(x,o)\,=\,\beta+\gamma\alpha,\quad\mathcal{B}(x,o)\,=\,-\gamma\alpha.$

Again, because the model is linear, $\mathcal{C}(x,o)$ is independent of
$p(U|x,o)$. If $\gamma\neq 0$ then $U_{V_{1}}\in\mathcal{U}^{C}$, and if
$\alpha\neq 0$ then $U_{V_{1}}\not\perp\\!\\!\\!\perp X|V_{1}$. Thus (2) is
violated and $\mathcal{B}(x,o)\neq 0$.

Selection bias. Given the model in Figure 4, we can compute

$\mathcal{A}(x,o)\,=\,\frac{\alpha-\gamma\beta}{1+\gamma^{2}},\quad\mathcal{C}(x,o)\,=\,\alpha,\quad\mathcal{B}(x,o)\,=\,-\frac{\gamma(\beta+\gamma\alpha)}{1+\gamma^{2}}.$

Once more, because the model is linear, $\mathcal{C}(x,o)$ is independent of
$p(U|x,o)$. If $\gamma\neq 0$ and either $\alpha\neq 0$ or $\beta\neq 0$, then
$U_{Y}\not\perp\\!\\!\\!\perp X|V_{1}$. Since $U_{Y}\in\mathcal{U}^{C}$, this
violates (2) and $\mathcal{B}(x,o)\neq 0$.

###  3.2 The lesser of two evils: covariate adjustment with missing data

In most real applications we cannot adjust for all confounders. For example,
take the causal graph in Figure 5. Whenever we model a mediator $V_{2}$, there
is likely to exist some variable $V_{1}$ for which we do not have data acting
as a confounder between $X$ and $V_{2}$. Theorem 2.4 says that in order to
achieve identifiability we would need to observe $V_{1}$ and not observe
$V_{2}$. However, as we cannot observe $V_{1}$, we are left with the following
question: Should we observe $V_{2}$ and incur overcontrol bias, or not observe
it and have confounding bias?

$V_{1}$$X$$V_{2}$$Y$$U_{V_{1}}$$U_{X}$$U_{V_{2}}$$U_{Y}$

$\displaystyle V_{1}$ $\displaystyle=U_{V_{1}}$ $\displaystyle X$
$\displaystyle=\alpha\exp(V_{1})+U_{X}$ $\displaystyle V_{2}$
$\displaystyle=\beta X+\gamma V_{1}^{2}+U_{V_{2}}$ $\displaystyle Y$
$\displaystyle=\delta V_{2}+U_{Y}$

Figure 5: Confounding vs. overcontrol bias. If $V_{1}$ cannot be observed,
should we observe $V_{2}$ or not?

Let us bring this question to the non-linear model in Figure 5, where we take
$\alpha,\beta,\gamma,\delta$ to be scalar parameters and
$U_{V_{1}},U_{X},U_{V_{2}},U_{Y}\sim\mathcal{N}(0,1)$. Here, the average
partial effect on the treated is given by $\mathcal{C}(x,o)=\beta\delta$
independently whether we observe $V_{2}$ or not, therefore it does not help us
towards a decision. Rather, we want to compute the absolute marginal causal
bias $|\mathcal{B}(x,o)|$ and pick the case where it is smaller. To this
purpose, here we first estimate $p(U_{\mathcal{L}}|x,o)$ via a Laplace
approximation as described in Section 2.3, using an ADAM optimizer (Kingma &
Ba, 2014) to compute the MAP; then we compose $p(U|x,o)$ as in (9).

Figure 6: Marginal distributions of $p(U|x,o)$ Figure 7: Comparisons of
$|\mathcal{B}(x,o)|$ over $x$, for arbitrary parameters
$\alpha,\beta,\gamma,\delta$

Figure 6 shows its marginal distributions for an arbitrary configuration of
parameters and for $x=1$, when $V_{2}$ is unobserved and when $V_{2}=2$.

Given approximate posterior samples, we can compute a Monte Carlo estimator of
the bias in (2.10). Figure 7 shows the estimated $|\mathcal{B}(x,o)|$ both
when $V_{2}$ is unobserved and when $V_{2}=2$, for two different
configurations of parameters $\alpha,\beta,\gamma,\delta$ and for
$x\in(-20,20)$. We can see that the estimated $|\mathcal{B}(x,o)|$ is constant
over $x$ when $V_{2}$ is observed, which is expected since for this model the
overcontrol bias equals $\beta\delta$. The confounding bias, however, is non-
linear in $x$. Just increasing $\alpha$ from $1$ to $5$, the range of values
of $x$ for which observing $V_{2}$ is better than not observing it changes
drastically. As a sanity check, we also make sure that the estimated value of
$\mathcal{B}(x,o)$ is close to $0$ if we could hypothetically observe only $X$
and $V_{1}$; this turns out in the order of $10^{-13}$, confirming that the
estimation procedure is sensible.

In conclusion, as expected the decision depends on specific observations and
parameters, potentially learned from data. If we can model missing
confounders, feasible bias estimation provides a quantitative tool to take
concrete decisions under missing data.

###  3.3 A simulated study of statins and atherosclerotic cardiovascular
disease

We follow the example in (Zivich & Breskin, 2021) and conduct a simulation
study of statins and subsequent atherosclerotic cardiovascular disease
(ASCVD). We adapt the data generation procedure as presented by the authors to
satisfy the requirements of Theorem 2.10. In particular, we assume a
continuous treatment variable $X$ representing the strength of each prescribed
medicine with support in the range $[0,1]$. This dose intensity of statin is
relative to the maximal dose, that is, 80 mg for atorvastatin and 40 mg for
rosuvastatin. This better reflects reality as a typical daily dose varies
depending on the patient characteristics and other confounders: age $A$, pre-
treatment low-density lipoprotein $L$, frailty $F$, diabetes $D$, and ASCVD
risk score $R$. We consider another measurement $M$ of low-density lipoprotein
that happens post-treatment. The outcome $Y$ indicates the observed incidence
of ASCVD. Finally, we record whether a patient had severe headache $H$, which
can be an adverse reaction to statin and/or caused by ASCVD.

We implement the data generation process as a $\theta$-parameterized joint
distribution of $X$ and $Y$ as well as covariates
$\mathcal{V}=\\{A,L,F,R,D,M,H\\}$. The model in use is complex and non-linear;
see Appendix C. Figure 8 shows some distributions and statistics recovered
from the data, which, although the model is not exactly the same, are
relatively close to those reported by the authors.

| size | age | diabetes | headache | log(pre-LDL) | log(post-LDL) | risk score | ASCVD
---|---|---|---|---|---|---|---|---
statin | 846 | 56.629 | 0.229 | 0.195 | 4.916 | 4.730 | 0.099 | 0.354
no statin | 2154 | 53.092 | 0.020 | 0.042 | 4.860 | 4.878 | 0.076 | 0.301

Figure 8: Distributions and statistics of the data generating process

| ${\mathcal{O}=\\{A\\!,\\!L\\!,\\!F\\!,\\!D\\}}$ | ${\mathcal{O}=\emptyset}$ | ${\mathcal{O}=\\{A\\!,\\!L\\!,\\!F\\!,\\!D\\!,\\!M\\}}$ | ${\mathcal{O}=\\{A\\!,\\!L\\!,\\!F\\!,\\!D\\!,\\!H\\}}$ | ${\mathcal{O}=\\{H\\!,\\!M\\}}$
---|---|---|---|---|---
effect | $-0.112\ \ \ (+0.045)$ | $-0.113\ \ \ (+0.004)$ | $-0.112\ \ \ (+0.045)$ | $-0.112\ \ \ (+0.047)$ | $-0.112\ \ \ (+0.027)$
bias | $+0.000\ \ \ (+0.000)$ | $+0.417\ \ \ (+0.895)$ | $-0.005\ \ \ (+0.005)$ | $-0.088\ \ \ (+0.048)$ | $+0.100\ \ \ (+0.608)$

Figure 9: Distributions and means (standard deviations) of average partial
effect on the treated and causal bias

In order to estimate $\mathcal{C}(x,o)$ and $\mathcal{B}(x,o)$, here we adopt
an importance sampling technique as described in Section 2.3, choosing the
prior $p(U_{\mathcal{L}})$ as importance density and $N=10^{5}$. Figure 9
shows distributions of estimated $\mathcal{C}(x,o)$ and $\mathcal{B}(x,o)$
with $x$ and $o$ randomly generated from the model, for different sets of
observed covariates $\mathcal{O}$; the table shows their mean estimates, with
standard deviation in brackets. When $\mathcal{O}=\\{A,L,F,D\\}$ the bias
distribution is peaked around zero. This is expected, since all confounders
are observed and only them, so the model should exhibit no bias. When
$\mathcal{O}=\emptyset$ we have confounding bias, which we see to be very
pronounced. In comparison to the case before, here the distribution of
$\mathcal{C}(x,o)$ is very peaked, showing the large impact of the
confounders. This is also the case studied in (Zivich & Breskin, 2021), where
a true causal effect of $-0.1082$ is reported in the authors’ code and is very
close to what we estimate for our modified model. When
$\mathcal{O}=\\{A,L,F,D,M\\}$ we have overcontrol bias, but this does not seem
particularly impactful. If $\mathcal{O}=\\{A,L,F,D,H\\}$ we have endogenous
selection bias, which is smaller than the confounding bias but stronger than
the overcontrol one. Finally, if $\mathcal{O}=\\{H,M\\}$ all biases interact
together. Interestingly, in this case the observation of post-treatment
effects such as $H$ and $M$, while themselves contributing to induce bias,
partially mitigates the effect of confounding bias.

##  4 Conclusion and future direction

In this work we developed a new complete criterion for identifiability of
causal effects via covariate adjustment. We provided a mathematical relation
between association, causal effect, and bias, showing the bias is zero if such
criterion holds. We developed a formulation of the bias where all terms can be
estimated from data using automatic differentiation. We demonstrated the
usefulness of our results over several experiments. A natural future work
direction is to exploit the bias characterization as a causal regularizer for
the training of highly non-linear machine learning models where we cannot
easily control for causal bias.

## References

  * Abrevaya et al. (2015) Abrevaya, J., Hsu, Y.-C., and Lieli, R. P.  Estimating conditional average treatment effects.  _Journal of Business & Economic Statistics_, 33(4):485–505, 2015.
  * Altshuler (1981) Altshuler, B.  Modeling of dose-response relationships.  _Environmental health perspectives_ , 42:23–27, 1981.
  * Angrist & Pischke (2008) Angrist, J. D. and Pischke, J.-S.  _Mostly harmless econometrics: An empiricist’s companion_.  Princeton university press, 2008.
  * Austin (2019) Austin, P. C.  Assessing covariate balance when using the generalized propensity score with quantitative or continuous exposures.  _Statistical methods in medical research_ , 28(5):1365–1377, 2019.
  * Bishop (2006) Bishop, C. M.  _Pattern Recognition and Machine Learning_.  Springer, 2006.  ISBN 978-0-387-31073-2.
  * Blei et al. (2017) Blei, D. M., Kucukelbir, A., and McAuliffe, J. D.  Variational inference: A review for statisticians.  _Journal of the American statistical Association_ , 112(518):859–877, 2017.
  * Brooks et al. (2011) Brooks, S., Gelman, A., Jones, G., and Meng, X.-L.  _Handbook of Markov Chain Monte Carlo_.  CRC press, 2011.
  * Chernozhukov et al. (2018) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J.  Double/debiased machine learning for treatment and structural parameters, 2018.
  * Chipman et al. (2010) Chipman, H. A., George, E. I., McCulloch, R. E., et al.  Bart: Bayesian additive regression trees.  _The Annals of Applied Statistics_ , 4(1):266–298, 2010.
  * Cinelli & Hazlett (2020) Cinelli, C. and Hazlett, C.  Making sense of sensitivity: Extending omitted variable bias.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 82(1):39–67, 2020.
  * Doucet et al. (2001) Doucet, A., De Freitas, N., and Gordon, N.  An introduction to sequential monte carlo methods.  In _Sequential Monte Carlo methods in practice_ , pp. 3–14. Springer, 2001.
  * Elwert (2013) Elwert, F.  Graphical causal models.  In _Handbook of causal analysis for social research_ , pp. 245–273. Springer, 2013.
  * Gain & Shpitser (2018) Gain, A. and Shpitser, I.  Structure learning under missing data.  In _International Conference on Probabilistic Graphical Models_ , pp. 121–132. PMLR, 2018.
  * Galagate (2016) Galagate, D.  _Causal inference with a continuous treatment and outcome: Alternative estimators for parametric dose-response functions with applications._ PhD thesis, 2016.
  * Gill & Robins (2001) Gill, R. D. and Robins, J. M.  Causal inference for complex longitudinal data: The continuous case.  _Annals of Statistics_ , pp. 1785–1811, 2001.
  * Hernán & Robins (2020) Hernán, M. A. and Robins, J. M.  Causal inference: What if, 2020.
  * Hill (2011) Hill, J. L.  Bayesian nonparametric modeling for causal inference.  _Journal of Computational and Graphical Statistics_ , 20(1):217–240, 2011.
  * Hines et al. (2021) Hines, O., Diaz-Ordaz, K., and Vansteelandt, S.  Parameterising the effect of a continuous exposure using average derivative effects, 2021.
  * Hirano & Imbens (2004) Hirano, K. and Imbens, G. W.  The propensity score with continuous treatments.  _Applied Bayesian modeling and causal inference from incomplete-data perspectives_ , 226164:73–84, 2004.
  * Imai & Van Dyk (2004) Imai, K. and Van Dyk, D. A.  Causal inference with general treatment regimes: Generalizing the propensity score.  _Journal of the American Statistical Association_ , 99(467):854–866, 2004.
  * Imbens (2000) Imbens, G. W.  The role of the propensity score in estimating dose-response functions.  _Biometrika_ , 87(3):706–710, 2000.
  * Imbens & Rubin (2015) Imbens, G. W. and Rubin, D. B.  _Causal inference in statistics, social, and biomedical sciences_.  Cambridge University Press, 2015.
  * Jordan et al. (1999) Jordan, M., Jaakkola, and Saul, L.  An Introduction to Variational Methods for Graphical Models.  _Machine Learning_ , (37):183–233, 1999.
  * Kennedy et al. (2017) Kennedy, E. H., Ma, Z., McHugh, M. D., and Small, D. S.  Nonparametric methods for doubly robust estimation of continuous treatment effects.  _Journal of the Royal Statistical Society. Series B, Statistical Methodology_ , 79(4):1229, 2017.
  * Kingma & Ba (2014) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Koller & Friedman (2009) Koller, D. and Friedman, N.  _Probabilistic Graphical Models: Principles and Techniques_.  Adaptive computation and machine learning. MIT Press, 2009.  ISBN 9780262013192.
  * Lauritzen (2001) Lauritzen, S. L.  Causal inference from graphical models.  _Complex stochastic systems_ , pp. 63–107, 2001.
  * Little & Rubin (2019) Little, R. J. and Rubin, D. B.  _Statistical analysis with missing data_ , volume 793.  John Wiley & Sons, 2019.
  * Murphy (2012) Murphy, K. P.  _Machine learning: A probabilistic perspective_.  MIT press, 2012.
  * Neal (2001) Neal, R. M.  Annealed importance sampling.  _Statistics and computing_ , 11(2):125–139, 2001\.
  * Papamakarios et al. (2019) Papamakarios, G., Nalisnick, E., Rezende, D. J., Mohamed, S., and Lakshminarayanan, B.  Normalizing flows for probabilistic modeling and inference.  _arXiv preprint arXiv:1912.02762_ , 2019.
  * Pearl (2009a) Pearl, J.  Causal inference in statistics: An overview.  _Statistics Surveys_ , 3:96–146, 2009a.  doi: 10.1214/09-SS057.
  * Pearl (2009b) Pearl, J.  _Causality_.  Cambridge university press, 2009b.
  * Pearl (2009c) Pearl, J.  _Causality: Models, Reasoning and Inference_.  Cambridge University Press, 2nd edition, 2009c.
  * Rezende & Mohamed (2015) Rezende, D. and Mohamed, S.  Variational inference with normalizing flows.  In _International Conference on Machine Learning_ , pp. 1530–1538. PMLR, 2015.
  * Rosenbaum & Rubin (1983) Rosenbaum, P. R. and Rubin, D. B.  The central role of the propensity score in observational studies for causal effects.  _Biometrika_ , 70(1):41–55, 1983.
  * Rosenbaum et al. (2010) Rosenbaum, P. R. et al.  _Design of observational studies_ , volume 10.  Springer, 2010.
  * Shpitser & Pearl (2006) Shpitser, I. and Pearl, J.  Identification of joint interventional distributions in recursive semi-Markovian causal models.  In _Proceedings of the 21st National Conference on Artificial Intelligence_ , volume 2, pp. 1219–1226. AAAI Press, 2006.
  * Shpitser et al. (2010) Shpitser, I., VanderWeele, T., and Robins, J. M.  On the validity of covariate adjustment for estimating causal effects.  In _Proceedings of the 26th Conference on Uncertainty in Artificial Intelligence, UAI 2010_ , pp. 527–536. AUAI Press, 2010.
  * Shun & McCullagh (1995) Shun, Z. and McCullagh, P.  Laplace approximation of high dimensional integrals.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 57(4):749–760, 1995.
  * Silva (2016) Silva, R.  Observational-interventional priors for dose-response learning.  In Lee, D., Sugiyama, M., Luxburg, U., Guyon, I., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 29, 2016\.
  * Spirtes et al. (2000) Spirtes, P., Glymour, C., and Scheines, R.  _Causation, Prediction, and Search_.  MIT press, 2nd edition, 2000.
  * Tran et al. (2016) Tran, D., Ruiz, F. J., Athey, S., and Blei, D. M.  Model criticism for Bayesian causal inference.  _arXiv preprint arXiv:1610.09037_ , 2016.
  * Wang (2015) Wang, J.  _Exposure-Response Modeling: Methods and Practical Implementation_ , volume 84.  CRC press, 2015.
  * Wu et al. (2018) Wu, X., Mealli, F., Kioumourtzoglou, M.-A., Dominici, F., and Braun, D.  Matching on generalized propensity scores with continuous exposures.  _arXiv preprint arXiv:1812.06575_ , 2018.
  * Zivich & Breskin (2021) Zivich, P. N. and Breskin, A.  Machine learning for causal inference: On the use of cross-fit estimators.  _Epidemiology_ , 32(3):393–401, 2021.

##  Appendix A Proofs

###  A.1 Proof of Theorem 2.4

Let us first report the adjustment criterion introduced in (Shpitser et al.,
2010), adapted to our framework where, for sake of simplicity, we required
$\mathcal{X}$ and $\mathcal{Y}$ to be singletons. We stress that, unlike the
definition in (Shpitser et al., 2010), by convention in this work we defined
parents, ancestors and descendants of a variable not to include the variable
itself (see (Lauritzen, 2001)).

######  Definition A.1.

Given a DAG $\mathcal{G}$, a set of nodes
$\mathcal{O}\in\mathcal{V}\setminus\\{X,Y\\}$ satisfies the adjustment
criterion relative to $(X,Y)$ in $\mathcal{G}$ if:

  * •

$\mathcal{O}$ does not include descendants of any variable laying on a causal
path from $X$ to $Y$, that is a path from $X$ to $Y$ where all arrows point
away from $X$;

  * •

$\mathcal{O}$ blocks all non-causal paths from $X$ to $Y$ in $\mathcal{G}$.

The adjustment criterion was shown in (Shpitser et al., 2010) to be complete
for identifiability of causal effects via covariate adjustment. We now show
that such criterion and (2) are equivalent. In the proof, we will use
graphical concepts such as blocked/unblocked path, fork and collider; see
(Pearl, 2009b) for their definitions.

Let us proceed by contradiction. Suppose there exists
$O_{j}\in\mathcal{O}\cap\text{De}(V_{i})\cup\\{V_{i}\\}$ with $V_{i}$ laying
on a causal path from $X$ to $Y$. Without loss of generality, suppose that a
causal path from $X$ to $V_{i}$ is unblocked, otherwise we could repeat the
argument for the first variable that blocks the path instead. Analogously,
without loss of generality suppose that a causal path form $V_{i}$ to $O_{j}$
is unblocked. Then we have $U_{V_{i}}\not\perp\\!\\!\\!\perp X|\mathcal{O}$.
Since $U_{V_{i}}\in\mathcal{U}^{C}$, this contradicts (2).

Now suppose instead there exists an unblocked non-causal path from $X$ to $Y$.
If the last arrow of such path points away from $Y$, then
$U_{Y}\not\perp\\!\\!\\!\perp X|\mathcal{O}$. Since $U_{Y}\in\mathcal{U}^{C}$,
this contradicts (2). Vice versa, if the last arrow of such path points into
$Y$, in order to be unblocked and non-causal the path must contain a fork
$V_{i}\not\in\mathcal{O}$ such that $V_{i}\in\text{An}(Y)$. Then
$U_{V_{i}}\not\perp\\!\\!\\!\perp X|\mathcal{O}$ and
$U_{V_{i}}\in\mathcal{U}^{C}$, contradicting (2).

This completes the proof that (2) implies the adjustment criterion. Let us now
show the other direction. By contradiction, suppose (2) does not hold. Then
there exists a $U_{V_{i}}\in\mathcal{U}\cap\text{An}^{+}(Y)$ such that
$U_{V_{i}}\not\perp\\!\\!\\!\perp Y|X,\mathcal{O}\setminus\text{De}(X)$ and
$U_{V_{i}}\not\perp\\!\\!\\!\perp X|\mathcal{O}$. Note that $V_{i}\neq X$
because $U_{X}\perp\\!\\!\\!\perp Y|X,\mathcal{O}\setminus\text{De}(X)$. Also,
we must have $V_{i}\in\text{An}(Y)$ or $V_{i}=Y$.

Suppose that $V_{i}=Y$. Then there exists some observed variable
$O_{j}\in\text{De}(Y)$ such that a path from $X$ to $O_{j}$ is unblocked,
otherwise $U_{Y}\perp\\!\\!\\!\perp X|\mathcal{O}$. Consequently, there exists
an unblocked path from $X$ to $Y$. If the latter is non-causal, this
contradicts the second item of the adjustment criterion. Otherwise, since also
the path from $Y$ to $O_{j}$ is causal by definition of descendant, we have a
contradiction with the first item instead.

Now suppose instead that $V_{i}\in\text{An}(Y)\setminus\mathcal{X}$ and
$V_{i}\in\mathcal{O}$. Since $U_{V_{i}}\not\perp\\!\\!\\!\perp X|\mathcal{O}$
there exists an unblocked path from $X$ to $V_{i}$ with last arrow pointing
into $V_{i}$. If, in addition, $V_{i}\in\mathcal{O}\setminus\text{De}(X)$,
then there exists an unblocked path from $V_{i}$ to $Y$ with first arrow into
$V_{i}$, otherwise $U_{V_{i}}\perp\\!\\!\\!\perp
Y|X,\mathcal{O}\setminus\text{De}(X)$. Then there is an unblocked non-causal
path from $X$ to $Y$ with observed collider $V_{i}$, which contradicts the
second item of the adjustment criterion. If, vice versa,
$V_{i}\in\mathcal{O}\cap\text{De}(X)$, then there is a causal path from $X$ to
$V_{i}$ and one from $V_{i}$ to $Y$ with $V_{i}$ observed, which contradicts
the first item instead.

On the other hand, suppose $V_{i}\in\text{An}(Y)\setminus\mathcal{X}$ and
$V_{i}\not\in\mathcal{O}$. Since $U_{V_{i}}\not\perp\\!\\!\\!\perp
X|\mathcal{O}$ there exists an unblocked path from $X$ to $V_{i}$ with last
arrow pointing away from $V_{i}$. If a causal path from $V_{i}$ to $Y$ is
unblocked, then there is an unblocked non-causal path from $X$ to $Y$ with
unobserved fork $V_{i}$, which contradicts the second item of the adjustment
criterion. If a causal path from $V_{i}$ to $Y$ is blocked by some variable in
$\mathcal{O}\cap\text{De}(X)$, then there is a blocked causal path from $X$ to
$Y$, which contradicts the first item. Otherwise, since
$U_{V_{i}}\in\mathcal{U}^{C}$ there must be a non-causal path from $V_{i}$ to
$Y$, with first arrow pointing away from $V_{i}$, that is unblocked when
observing variables in $\mathcal{O}\setminus\text{De}(X)$. If this path is
also unblocked by $\mathcal{O}$, then there is a contradiction with the second
item of the adjustment criterion. Otherwise, take the last variable
$O_{j}\in\mathcal{O}\cap\text{De}(X)$ that is blocking the path. This can only
happen if $O_{j}\in\text{An}(Y)$, because otherwise there would need to be
$O_{k}\in\mathcal{O}\cap\text{De}(O_{j})$ that is either an observed collider
or an observed descendant of a collider, which would make the path from
$V_{i}$ to $Y$ blocked when observing variables in
$\mathcal{O}\setminus\text{De}(X)$. Hence we have a contradiction with the
first item of the adjustment criterion and conclude the proof.

###  A.2 Proof of Proposition 2.7

Conditioned over $X=x$ and $O=o$, we have $Y=f_{Y}=f_{Y}^{x,\underline{o}}$,
where $=$ stands for equality in distribution. By Definition 2.1, the causal
exogenous set $\mathcal{U}^{C}$ uniquely determines $f_{Y}^{x,\underline{o}}$,
as it includes all exogenous random variables that are ancestors of $Y$ and
are not independent of $Y$ given $x$ and $\underline{o}$. Then by the Law of
the Unconscious Statistician we can write

$\mathbb{E}_{Y|x,o}[Y]=\int y\,p(y|x,o)\,dy=\int
f_{Y}^{x,\underline{o}}\,p(u^{c}|x,o)\,du^{c}.$

Next, as we assumed $f_{Y}$ differentiable with respect to $X$ and satisfying
sufficient regularity assumptions, by integral Leibniz rule we have

$\displaystyle\nabla_{x}\mathbb{E}_{Y|x,o}[Y]$
$\displaystyle=\int\nabla_{x}\left(f_{Y}^{x,\underline{o}}\,p(u^{c}|x,o)\right)\,du^{c}$
$\displaystyle=\int\left(\nabla_{x}f_{Y}^{x,\underline{o}}+f_{Y}^{x,\underline{o}}\,\nabla_{x}\log
p(u^{c}|x,o)\right)p(u^{c}|x,o)\,du^{c}$
$\displaystyle=\mathbb{E}_{U^{C}|x,o}[\nabla_{x}\,f_{Y}^{x,\underline{o}}]+\mathbb{E}_{U^{C}|x,o}[f_{Y}^{x,\underline{o}}\,\nabla_{x}\log
p(U^{C}|x,o)],$

where we used product rule and the fact that
$\nabla_{x}p(u^{c}|x,o)=p(u^{c}|x,o)\nabla_{x}\log p(u^{c}|x,o)$. Finally, if
the causal effect of $X$ on $Y$ can be identified via covariate adjustment
then Corollary 2.5 holds, which trivially implies $\mathcal{B}(x,o)=0$.

###  A.3 Proof of Theorem 2.9

First, we observe that $\nabla_{x}\log p(U^{C}|x,o)=\nabla_{x}\log
p(U^{C},x,o)-\nabla_{x}\log p(x,o)$. Then we have

$\displaystyle\nabla_{x}\log p(x,o)$
$\displaystyle=\frac{1}{p(x,o)}\nabla_{x}p(x,o)$
$\displaystyle=\frac{1}{p(x,o)}\nabla_{x}\int p(u^{c},x,o)\,dc$
$\displaystyle=\frac{1}{p(x,o)}\int\nabla_{x}p(u^{c},x,o)\,dc$
$\displaystyle=\frac{1}{p(x,o)}\int\nabla_{x}\log
p(u^{c},x,o)p(u^{c},x,o)\,dc$ $\displaystyle=\int\nabla_{x}\log
p(u^{c},x,o)p(u^{c}|x,o)\,dc$
$\displaystyle=\mathbb{E}_{U^{C}|x,o}[\nabla_{x}\log p(U^{C},x,o)].$

Thus the bias can be rewritten as

$\displaystyle\mathcal{B}(x,o)$
$\displaystyle=\mathbb{E}_{U^{C}|x,o}[f_{Y}^{x,\underline{o}}\,\nabla_{x}\log
p(U^{C}|x,o)]$
$\displaystyle=\mathbb{E}_{U^{C}|x,o}[f_{Y}^{x,\underline{o}}\,\nabla_{x}\log
p(U^{C},x,o)]-\mathbb{E}_{U^{C}|x,o}[f_{Y}^{x,\underline{o}}\,\nabla_{x}\log
p(x,o)]$
$\displaystyle=\mathbb{E}_{U^{C}|x,o}[f_{Y}^{x,\underline{o}}\,\nabla_{x}\log
p(U^{C},x,o)]-\mathbb{E}_{U^{C}|x,o}[f_{Y}^{x,\underline{o}}]\,\nabla_{x}\log
p(x,o)$
$\displaystyle=\mathbb{E}_{U^{C}|x,o}[f_{Y}^{x,\underline{o}}\,\nabla_{x}\log
p(U^{C},x,o)]-\mathbb{E}_{U^{C}|x,o}[f_{Y}^{x,\underline{o}}]\,\mathbb{E}_{C|x,o}[\nabla_{x}\log
p(U^{C},x,o)]$
$\displaystyle=\text{Cov}_{U^{C}|x,o}\left(f_{Y}^{x,\underline{o}},\nabla_{x}\log
p(U^{C},x,o)\right).$

###  A.4 Proof of Theorem 2.10

Let us start by the expression of the bias in (7). We first note that
$\bar{f}_{Y}^{x,\underline{o}}=\int\bar{f}_{Y}^{x}p(u_{\underline{O}}|u^{c},\underline{o})\,du_{\underline{O}}$,
where $p(u_{\underline{O}}|u^{c},\underline{o})$ is a Dirac delta
distribution. In particular the latter does not depend on $x$ by definition of
$\underline{\c{O}}$. By further marginalization, this also implies that
$\mathbb{E}_{U^{C}|x,o}[f_{Y}^{x,\underline{o}}]=\mathbb{E}_{U|x,o}[f_{Y}^{x}]$.
For sake of notation, we will denote
$\bar{f}_{Y}^{x,\underline{o}}=f_{Y}^{x,\underline{o}}-\mathbb{E}_{U^{C}|x,o}[f_{Y}^{x,\underline{o}}]$
and $\bar{f}_{Y}^{x}=f_{Y}^{x}-\mathbb{E}_{U|x,o}[f_{Y}^{x}]$. Also, denote
$u=[u^{c},u_{\underline{O}},\tilde{u}]$, where $u^{c}$ is a realization of the
causal exogenous set and $\tilde{u}$ contains all variables in $u$ that are
neither in $u^{c}$ nor in $u_{\underline{O}}$. We have

$\displaystyle\mathcal{B}(x,o)$
$\displaystyle=\text{Cov}_{U^{C}|x,o}(f_{Y}^{x,\underline{o}},\nabla_{x}\log
p(u^{c},x,o))$
$\displaystyle=\mathbb{E}_{U^{C}|x,o}[\bar{f}_{Y}^{x,\underline{o}}\,\nabla_{x}\log
p(u^{c},x,o)]$
$\displaystyle=\frac{1}{p(x,o)}\int\bar{f}_{Y}^{x,\underline{o}}\,\nabla_{x}p(u^{c},x,o)\,du6c$
$\displaystyle=\frac{1}{p(x,o)}\int\left(\int\bar{f}_{Y}^{x}\,p(u_{\underline{O}}|u^{c},\underline{o})\,du_{\underline{O}}\right)\nabla_{x}\left(\int
p(u^{c},\tilde{u},x,o)\,d\tilde{u}\right)\,du^{c}$
$\displaystyle=\frac{1}{p(x,o)}\int\bar{f}_{Y}^{x}\nabla_{x}(p(u_{\underline{O}}|u^{c},\underline{o})\,p(u^{c},\tilde{u},x,o))\,du$
$\displaystyle=\frac{1}{p(x,o)}\int\bar{f}_{Y}^{x}\nabla_{x}p(u,x,o)\,du$

where (i) we marginalized $p(u^{c},x,o)$ over variables $\tilde{u}$; (ii) we
moved the integrals outside, which is possible because $\bar{f}_{Y}^{x}$ does
not depend on $\tilde{u}$ by definition of $\mathcal{U}^{C}$; (iii) we moved
$p(u_{\underline{O}}|u^{c},\underline{o})$ inside the gradient with respect to
$x$, which is allowed since it is independent of it; (iv) we rewrote the
density factorization as $p(u,x,o)$, that is correct because
$p(u_{\underline{O}}|u^{c},\underline{o})=p(u_{\underline{O}}|u^{c},x,o,\tilde{u})$
by definition of $\underline{\mathcal{O}}$.

Again for sake of notation, let us denote $z=[x,o]$. We indicate by $z_{i}$ a
generic variable $x$ or $o_{i}$, whereas $z_{-i}$ corresponds to all variables
in $z$ except $z_{i}$. Without loss of generality, we further order the
indices $i$ from parents to children, so that a variable $z_{i}$ may depend on
$z_{i-1}$ but not on $z_{i+1}$. We write $z_{<i}=[z_{1},\dots,z_{i-1}]$. Then
we can factorize

$p(u,z)=p(u)\prod_{i}p(z_{i}|z_{<i},u)=p(u)\prod_{i}\delta(g_{z_{i}}),$

where $\delta$ denotes a Dirac delta distribution and we define
$g_{z_{i}}=f_{Z_{i}}^{x,o}-z_{i}$. Replacing in the derivation above, we have

$\displaystyle\frac{1}{p(x,o)}\int\bar{f}_{Y}^{x}\,\nabla_{x}p(u,x,o)\,du$
$\displaystyle=\frac{1}{p(x,o)}\sum_{i}\int\bar{f}_{Y}^{x}\,p(u)p(z_{-i}|u)\nabla_{x}\delta(g_{z_{i}})\,du$
$\displaystyle=\frac{1}{p(x,o)}\sum_{i}\int\bar{f}_{Y}^{x}\,p(u)p(z_{-i}|u)\nabla_{g_{z_{i}}}\delta(g_{z_{i}})\nabla_{x}g_{z_{i}}\,du,$

where in the last step we used chain rule. Let us now focus on each term of
the sum. We decompose $u=[u_{Z_{i}},u_{-{Z_{i}}}]$, where $u_{-{Z_{i}}}$
contains all variables in $u$ but $u_{Z_{i}}$. Observe that
$p(z_{-i}|u)=p(z_{-i}|u_{-Z_{i}})$ because $z_{-i}$ is independent of
$u_{Z_{i}}$. We then proceed with a change of variable, i.e.
$p(u_{Z_{i}})\,du_{Z_{i}}=p(g_{z_{i}})\,dg_{z_{i}}$. Replacing above, we get

$\displaystyle\int\bar{f}_{Y}^{x}\,p(u)p(z_{-i}|u)\nabla_{g_{z_{i}}}\delta(g_{z_{i}})\nabla_{x}g_{z_{i}}\,du$
$\displaystyle=\int\bar{f}_{Y}^{x}\,p(u_{Z_{i}})p(u_{-Z_{i}})p(z_{-i}|u_{-Z_{i}})\nabla_{g_{z_{i}}}\delta(g_{z_{i}})\nabla_{x}g_{z_{i}}\,du_{Z_{i}}du_{-Z_{i}}$
$\displaystyle=\int\bar{f}_{Y}^{x}\,p(g_{z_{i}})p(u_{-Z_{i}})p(z_{-i}|u_{-Z_{i}})\nabla_{g_{z_{i}}}\delta(g_{z_{i}})\nabla_{x}g_{z_{i}}\,dg_{z_{i}}du_{-Z_{i}}$
$\displaystyle=-\int\nabla_{g_{z_{i}}}\left(\bar{f}_{Y}^{x}\,p(g_{z_{i}})\nabla_{x}g_{z_{i}}p(u_{-Z_{i}})p(z_{-i}|u_{-Z_{i}})\right)\delta(g_{z_{i}})\,dg_{z_{i}}du_{-Z_{i}},$

where in the last step we used the definition of derivative of a Dirac delta,
i.e. $\int h(a)\nabla_{a}\delta(a)\,da=-\int\nabla_{a}h(a)\delta(a)\,da$,
where $h$ is a smooth function. This is possible because we assumed $f_{Y}$ to
be differentiable in $X$ and $V_{i}\in\mathcal{O}$. Note that
$\nabla_{g_{z_{i}}}\nabla_{x}g_{z_{i}}=\nabla_{x}\nabla_{g_{z_{i}}}g_{z_{i}}=0$;
furthermore, $\nabla_{g_{z_{i}}}p(u_{-Z_{i}})=0$ and
$\nabla_{g_{z_{i}}}p(z_{-i}|u_{-Z_{i}})=0$ because they do not depend on
$u_{Z_{i}}$. Also, notice that

$\nabla_{g_{z_{i}}}p(g_{z_{i}})=p(g_{z_{i}})\nabla_{g_{z_{i}}}\log
p(g_{z_{i}})=p(g_{z_{i}})\left(\nabla_{g_{z_{i}}}\log
p(u_{Z_{i}})-\nabla_{g_{z_{i}}}\log|\det(\nabla_{u_{Z_{i}}}g_{z_{i}})|\right).$

We can always express
$\det(\nabla_{u_{Z_{i}}}g_{z_{i}})=\sum_{\rho}\text{sign}(\rho)\prod_{k}\partial_{u_{Z_{i}}}[g_{z_{i}}]_{\rho(k)}$,
where $\rho$ is a permutation over the indices,
$\text{sign}(\rho)\in\\{+1,-1\\}$ is its sign and $[g_{z_{i}}]_{\rho(k)}$
denotes the ${\rho(k)}$-component of $g_{z_{i}}$. Because
$\partial_{u_{Z_{i}}}\nabla_{g_{z_{i}}}[g_{z_{i}}]_{\rho(i)}=0$, we conclude
that $\nabla_{g_{z_{i}}}\log|\det(\nabla_{u_{Z_{i}}}g_{z_{i}})|)=0$ and
$\nabla_{g_{z_{i}}}p(g_{z_{i}})=p(g_{z_{i}})\nabla_{g_{z_{i}}}\log
p(u_{Z_{i}})$. By replacing in the derivation above, we have

$\displaystyle-\int\nabla_{g_{z_{i}}}\left(\bar{f}_{Y}^{x}\,p(g_{z_{i}})\nabla_{x}g_{z_{i}}\,p(u_{-Z_{i}})p(z_{-i}|u_{-Z_{i}})\right)\delta(g_{z_{i}})\,dg_{z_{i}}du_{-Z_{i}}$
$\displaystyle=-\int\left(\nabla_{g_{z_{i}}}\bar{f}_{Y}^{x}+\bar{f}_{Y}^{x}\,\nabla_{g_{z_{i}}}\log
p(u_{Z_{i}})\right)\nabla_{x}g_{z_{i}}\,p(g_{z_{i}})p(u_{-Z_{i}})p(z_{-i}|u_{-Z_{i}})\delta(g_{z_{i}})\,dg_{z_{i}}du_{-Z_{i}}$
$\displaystyle=-\int\left(\nabla_{g_{z_{i}}}\bar{f}_{Y}^{x}+\bar{f}_{Y}^{x}\,\nabla_{g_{z_{i}}}\log
p(u_{Z_{i}})\right)\nabla_{x}g_{z_{i}}\,p(u_{Z_{i}})p(u_{-Z_{i}})p(z_{-i}|u_{-Z_{i}})\delta(g_{z_{i}})\,du$
$\displaystyle=-\int\left(\nabla_{g_{z_{i}}}\bar{f}_{Y}^{x}+\bar{f}_{Y}^{x}\,\nabla_{g_{z_{i}}}\log
p(u_{Z_{i}})\right)\nabla_{x}g_{z_{i}}\,p(u,z)\,du,$

where in the last two steps we undid the change of variables and reassembled
the joint density. We now sum together all terms over $i$, scale by $p(z)$ and
remember that we defined $z=[x,o]$,
$\bar{f}_{Y}^{x}=f_{Y}^{x}-\mathbb{E}_{U|x,o}[f_{Y}^{x}]$ and
$g_{v_{i}}=f_{V_{i}}^{x,o}-v_{i}$ for $V_{i}\in\mathcal{X}\cup\mathcal{O}$. We
have

$\displaystyle\mathcal{B}(x,o)$
$\displaystyle=-\frac{1}{p(z)}\sum_{i}\int\left(\nabla_{g_{z_{i}}}\bar{f}_{Y}^{x}+\bar{f}_{Y}^{x}\,\nabla_{g_{z_{i}}}\log
p(u_{Z_{i}})\right)\nabla_{x}g_{z_{i}}\,p(u,z)\,du$
$\displaystyle=-\sum_{i}\int\left(\nabla_{g_{z_{i}}}\bar{f}_{Y}^{x}+\bar{f}_{Y}^{x}\,\nabla_{g_{z_{i}}}\log
p(u_{Z_{i}})\right)\nabla_{x}g_{z_{i}}\,p(u|z)\,du$
$\displaystyle=-\sum_{V_{i}\in X\cup
O}\int\left(\nabla_{g_{v_{i}}}\bar{f}_{Y}^{x}+\bar{f}_{Y}^{x}\,\nabla_{g_{v_{i}}}\log
p(u_{V_{i}})\right)\nabla_{x}g_{v_{i}}\,p(u|x,o)\,du$
$\displaystyle=-\sum_{V_{i}\in X\cup
O}\mathbb{E}_{U|x,o}\left[\left(\nabla_{f_{V_{i}}^{x,o}}f_{Y}^{x}+\left(f_{Y}^{x}-\mathbb{E}_{U|x,o}[f_{Y}^{x}]\right)\nabla_{f_{V_{i}}^{x,o}}\log
p(u_{V_{i}})\right)\nabla_{x}(f_{V_{i}}^{x,o}-v_{i})\right],$

where we named $V_{i}$ to be a generic variable in
$\mathcal{X}\cup\mathcal{O}$ and we expressed the integral as an expectation
under the probability density $p(u|x,o)$. We can finally observe that for a
generic function $h$ differentiable in $f_{V_{i}}^{x,o}$, by chain rule we
have
$\nabla_{u_{V_{i}}}h=\nabla_{f_{V_{i}}^{x,o}}h\nabla_{u_{V_{i}}}f_{V_{i}}^{x,o}$.
Because $f_{V_{i}}^{x,o}$ is invertible with respect to $U_{V_{i}}$ for
$V_{i}\in\mathcal{X}\cup\mathcal{O}$, the Jacobian
$\nabla_{u_{V_{i}}}f_{V_{i}}^{x,o}$ is also invertible and we can write
$\nabla_{f_{V_{i}}^{x,o}}h=\nabla_{u_{V_{i}}}h\left(\nabla_{u_{V_{i}}}f_{V_{i}}^{x,o}\right)^{-1}$.
Replacing in the derivation above, we conclude that

$\mathcal{B}(x,o)=-\\!\\!\\!\sum_{V_{i}\in\mathcal{X}\cup\mathcal{O}}\\!\\!\mathbb{E}_{U|x,o}\Big{[}\Big{(}\nabla_{u_{V_{i}}}f_{Y}^{x}+\left(f_{Y}^{x}-\mathbb{E}_{U|x,o}[f_{Y}^{x}]\right)\nabla_{u_{V_{i}}}\\!\log
p(U_{V_{i}})\Big{)}(\nabla_{u_{V_{i}}}f_{V_{i}}^{x,o})^{-1}\,\nabla_{x}(f_{V_{i}}^{x,o}-v_{i})\Big{]}.$

##  Appendix B Confounding, overcontrol and endogenous selection bias

###  B.1 Confounding bias

Let us consider the model in Figure 4. Following (6), the average partial
effect on the treated can be expressed as

$\mathcal{C}(x,o)=\mathbb{E}_{U|x,o}[\nabla_{x}f_{Y}^{x}]=\mathbb{E}_{U|x,o}[\beta]=\beta.$

For the marginal causal bias, we have $\nabla_{u_{X}}f_{Y}^{x,o}=0$,
$\nabla_{u_{X}}\log p(U_{X})=-U_{X}$, $\nabla_{u_{X}}f_{X}^{x,o}=1$ and
$\nabla_{x}(f_{X}^{x,o}-x)=-1$. Then (2.10) gives

$\mathcal{B}(x,o)=-\mathbb{E}_{U|x,o}[\left(f_{Y}^{x}-\mathbb{E}_{U|x,o}[f_{Y}^{x}]\right)U_{X}].$

Notice that $p(u_{V_{1}},u_{X},u_{Y}|x)=p(u_{V_{1}},u_{X}|x)p(u_{Y})$, which
implies that $U_{Y}|X=x$ and $U_{X}|X=x$ are independent. This yields

$\displaystyle-\mathbb{E}_{U|x,o}[\left(f_{Y}^{x}-\mathbb{E}_{U|x,o}[f_{Y}^{x}]\right)U_{X}]$
$\displaystyle=-\mathbb{E}_{U|x,o}[(\beta x+U_{Y}-\mathbb{E}_{U|x,o}[\beta
x+U_{Y}])U_{X}]-\mathbb{E}_{U|x,o}[(\gamma U_{V_{1}}-\mathbb{E}_{U|x,o}[\gamma
U_{V_{1}}])U_{X}]$ $\displaystyle=-\mathbb{E}_{U|x,o}[\beta
x+U_{Y}-\mathbb{E}_{U|x,o}[\beta
x+U_{Y}]]\,\mathbb{E}_{U|x,o}[U_{X}]-\gamma\,\mathbb{E}_{U|x,o}[(U_{V_{1}}-\mathbb{E}_{U|x,o}[U_{V_{1}}])U_{X}]$
$\displaystyle=-\gamma\,\mathbb{E}_{U|x,o}[(U_{V_{1}}-\mathbb{E}_{U|x,o}[U_{V_{1}}])U_{X}]$
$\displaystyle=-\gamma\,\text{Cov}_{U|x,o}(U_{V_{1}},U_{X}),$

since $\mathbb{E}_{U|x,o}[\beta x+U_{Y}-\mathbb{E}_{U|x,o}[\beta x+U_{Y}]]=0$.
With some simple calculation, one can show that
$p(u_{X}|x)=N(\tfrac{x}{1+\alpha^{2}},\frac{\alpha^{2}}{1+\alpha^{2}})(u_{X})$
and $p(u_{V_{1}}|x)=N(\tfrac{\alpha
x}{1+\alpha^{2}},\tfrac{1}{1+\alpha^{2}})(u_{V_{1}})$, where
$N(\mu,\sigma^{2})(a)$ generally denotes a Gaussian probability density with
mean $\mu$, variance $\sigma^{2}$, evaluated at $a$. We then have that
$\mathbb{E}_{U|x,o}[U_{X}]=\tfrac{x}{1+\alpha^{2}}$,
$\mathbb{E}_{U|x,o}[U_{V_{1}}]=\tfrac{\alpha x}{1+\alpha^{2}}$ and
$\text{Var}_{U|x,o}(U_{V_{1}})=\tfrac{1}{1+\alpha^{2}}$. Furthermore, we have

$\displaystyle\mathbb{E}_{U|x,o}[U_{V_{1}}U_{X}]$ $\displaystyle=\int
u_{V_{1}}u_{X}p(u_{V_{1}},u_{X}|x)\,du_{X}\,du_{V_{1}}$ $\displaystyle=\int
u_{V_{1}}\int u_{X}p(u_{X}|u_{V_{1}},x)\,du_{X}\,p(u_{V_{1}}|x)\,du_{V_{1}}$
$\displaystyle=\int u_{V_{1}}\int u_{X}\delta(u_{X}+\alpha
u_{V_{1}}-x)\,du_{X}\,p(u_{V_{1}}|x)\,du_{V_{1}}$ $\displaystyle=\int
u_{V_{1}}(x-\alpha u_{V_{1}})\,p(u_{V_{1}}|x)du_{V_{1}}$
$\displaystyle=x\mathbb{E}_{U|x,o}[U_{V_{1}}]-\alpha\mathbb{E}_{U|x,o}[U_{V_{1}}^{2}]$
$\displaystyle=x\mathbb{E}_{U|x,o}[U_{V_{1}}]-\alpha(\text{Var}_{U|x,o}(U_{V_{1}})+\mathbb{E}_{U|x,o}[U_{V_{1}}]^{2}).$

Then we have all the elements to compute

$\mathcal{B}(x,o)=-\gamma\left(\mathbb{E}_{U|x,o}[U_{V_{1}}U_{X}]-\mathbb{E}_{U|x,o}[U_{V_{1}}]\mathbb{E}_{U|x,o}[U_{X}]\right)=\frac{\gamma\alpha}{1+\alpha^{2}}.$

Thus

$\mathcal{C}(x,o)+\mathcal{B}(x,o)=\beta+\frac{\gamma\alpha}{1+\alpha^{2}}.$

Let us check that the latter matches the association $\mathcal{A}(x,o)$ when
computed directly. We have

$\displaystyle\mathbb{E}_{Y|x,o}[Y]$ $\displaystyle=\beta
x+\gamma\mathbb{E}_{V_{1}|x,o}[V_{1}]+\mathbb{E}_{U|x,o}[U_{Y}]$
$\displaystyle=\beta x+\gamma\mathbb{E}_{U|x,o}[U_{V_{1}}]$
$\displaystyle=\beta x+\frac{\gamma\alpha}{1+\alpha^{2}}x.$

Then
$\mathcal{A}(x,o)=\nabla_{x}\mathbb{E}_{Y|x,o}[Y]=\beta+\tfrac{\gamma\alpha}{1+\alpha^{2}}$,
which verifies the statement.

###  B.2 Overcontrol bias

Let us consider the model in Figure 4. Following (6), we have

$\mathcal{C}(x,o)=\mathbb{E}_{U|x,o}[\beta+\gamma\alpha]=\beta+\gamma\alpha.$

For the bias, we have $\nabla_{u_{X}}f_{Y}^{x,o}=0$, $\nabla_{u_{X}}\log
p(U_{X})=-U_{X}$, $\nabla_{u_{X}}f_{X}^{x,o}=1$ and
$\nabla_{x}(f_{X}^{x,o}-x)=-1$. Then the contribution given by $X$ to the bias
is null, since

$\mathbb{E}_{U|x,o}[\left(f_{Y}^{x}-\mathbb{E}_{U|x,o}[f_{Y}^{x}]\right)U_{X}]=\mathbb{E}_{U|x,o}[f_{Y}^{x}-\mathbb{E}_{U|x,o}[f_{Y}^{x}]]\mathbb{E}_{U|x,o}[U_{X}]=0,$

where we used that $f_{Y}^{x}$ is independent of $U_{X}$ given $X=x$.

Regarding the contribution to the bias given by $V_{1}$, we have
$\nabla_{u_{V_{1}}}f_{Y}^{x}=\gamma$, $\nabla_{x}f_{V_{1}}^{x,o}=\alpha$ and
$\nabla_{u_{V_{1}}}f_{V_{1}}^{x,o}=1$. Moreover, $U_{Y}|V_{1}$ is independent
of $X$. Then

$\displaystyle\mathcal{B}(x,o)$
$\displaystyle=-\gamma\alpha+\alpha\mathbb{E}_{U|x,o}[\left(f_{Y}^{x}-\mathbb{E}_{U|x,o}[f_{Y}^{x}]\right)U_{V_{1}}]$
$\displaystyle=-\gamma\alpha+\alpha\mathbb{E}_{U|x,o}[\beta x+\gamma
v_{1}+U_{Y}-\mathbb{E}_{U|x,o}[\beta x+\gamma
v_{1}+U_{Y}]]\mathbb{E}_{U|x,o}[U_{V_{1}}]$ $\displaystyle=-\gamma\alpha$

since $\mathbb{E}_{U|x,o}[\beta x+\gamma v_{1}+U_{Y}-\mathbb{E}_{U|x,o}[\beta
x+\gamma v_{1}+U_{Y}]]=0$. Then

$\mathcal{C}(x,o)+\mathcal{B}(x,o)=\beta+\gamma\alpha-\gamma\alpha=\beta.$

Let us check the result by computing the marginal association directly. We
have $p(y|x,v_{1})=N(\beta x+\gamma v_{1},1)(y)$, whence
$\mathbb{E}_{Y|x,o}[Y]=\beta x+\gamma v_{1}$ and
$\mathcal{A}(x,o)=\nabla_{x}\mathbb{E}_{Y|x,o}[Y]=\beta$.

###  B.3 Endogenous selection bias

Let us consider the model in Figure 4. Following (6), we have

$\mathcal{C}(x,o)=\mathbb{E}_{U|x,o}[\nabla_{x}f_{Y}^{x}]=\mathbb{E}_{U|x,o}[\alpha]=\alpha.$

For the bias, we can argue that the contribution of $X$ is null for exactly
the same reason as in the overcontrol bias case (see Appendix B.2). Regarding
the contribution of $V_{1}$, we have $\nabla_{u_{V_{1}}}f_{Y}^{x}=0$,
$\nabla_{u_{V_{1}}}\log p(U_{V_{1}})=-U_{V_{1}}$,
$\nabla_{u_{V_{1}}}f_{V_{1}}^{x,o}=1$ and
$\nabla_{x}(f_{V_{1}}^{x,o}-{v_{1}})=\beta+\gamma\alpha$. Then

$\displaystyle\mathcal{B}(x,o)$
$\displaystyle=(\beta+\gamma\alpha)\mathbb{E}_{U|x,o}[\left(f_{Y}^{x}-\mathbb{E}_{U|x,o}[f_{Y}^{x}]\right)U_{V_{1}}]$
$\displaystyle=(\beta+\gamma\alpha)\mathbb{E}_{U|x,o}[\left(U_{Y}-\mathbb{E}_{U|x,o}[U_{Y}]\right)U_{V_{1}}]$
$\displaystyle=(\beta+\gamma\alpha)\text{Cov}_{U|x,o}(U_{Y},U_{V_{1}}).$

One can compute
$p(u_{Y}|x,v_{1})=N(\tfrac{\gamma}{1+\gamma^{2}}(v_{1}-(\beta+\gamma\alpha)x),\tfrac{1}{1+\gamma^{2}})(u_{Y})$
and
$p(u_{V_{1}}|x,v_{1})=N(\tfrac{1}{1+\gamma^{2}}(v_{1}-(\beta+\gamma\alpha)x),\tfrac{\gamma^{2}}{1+\gamma^{2}})(u_{V_{1}})$,
which respectively imply
$\mathbb{E}_{U|x,o}[U_{Y}]=\tfrac{\gamma}{1+\gamma^{2}}(v_{1}-\beta-\gamma\alpha)x$,
$\mathbb{E}_{U|x,o}[U_{V_{1}}]=\tfrac{x}{1+\gamma^{2}}(v_{1}-\beta-\gamma\alpha)$
and $\text{Var}_{U|x,o}(U_{V_{1}})=\tfrac{\gamma^{2}}{1+\gamma^{2}}$. If
$\gamma=0$, then $\mathbb{E}_{U|x,o}[U_{V_{1}}U_{Y}]=0$. Otherwise

$\displaystyle\mathbb{E}_{U|x,o}[U_{V_{1}}U_{Y}]$ $\displaystyle=\int
u_{V_{1}}u_{Y}p(u_{V_{1}},u_{Y}|x,v_{1})\,du_{Y}\,du_{V_{1}}$
$\displaystyle=\int u_{V_{1}}\int
u_{Y}p(u_{Y}|u_{V_{1}},x,v_{1})\,du_{Y}\,p(u_{V_{1}}|x,v_{1})\,du_{V_{1}}$
$\displaystyle=\int u_{V_{1}}\int u_{Y}\delta((\beta+\gamma\alpha)x+\gamma
u_{Y}+u_{V_{1}}-v_{1})\,du_{Y}\,p(u_{V_{1}}|x,v_{1})\,du_{V_{1}}$
$\displaystyle=\int
u_{V_{1}}\left(\frac{1}{\gamma}(v_{1}-(\beta+\gamma\alpha)x)-\frac{1}{\gamma}u_{V_{1}}\right)\,du_{V_{1}}$
$\displaystyle=\frac{1}{\gamma}(v_{1}-(\beta+\gamma\alpha)x)\mathbb{E}_{U|x,o}[U_{V_{1}}]-\frac{1}{\gamma}\mathbb{E}_{U|x,o}[U_{V_{1}}^{2}]$
$\displaystyle=\frac{1}{\gamma}(v_{1}-(\beta+\gamma\alpha)x)\mathbb{E}_{U|x,o}[U_{V_{1}}]-\frac{1}{\gamma}(\text{Var}_{U|x,o}(U_{V_{1}})+\mathbb{E}_{U|x,o}[U_{V_{1}}]^{2}).$

Then we have all the elements to compute

$\mathcal{B}(x,o)=(\beta+\gamma\alpha)\left(\mathbb{E}_{U|x,o}[U_{V_{1}}U_{Y}]-\mathbb{E}_{U|x,o}[U_{V_{1}}]\mathbb{E}_{U|x,o}[U_{Y}]\right)=-\frac{\gamma(\beta+\gamma\alpha)}{1+\gamma^{2}}.$

Thus

$\mathcal{C}(x,o)+\mathcal{B}(x,o)=\frac{\alpha-\gamma\beta}{1+\gamma^{2}}.$

Once more, we check that this is the same as computing the association
directly. Since $p(y|x,v_{1})=N(\tfrac{1}{1+\gamma^{2}}\left(\gamma
v_{1}+(\alpha-\beta\gamma)x\right),\tfrac{1}{1+\gamma^{2}})(y)$, we have
$\mathbb{E}_{Y|x,o}[Y]=\tfrac{1}{1+\gamma^{2}}\left(\gamma
v_{1}+(\alpha-\beta\gamma)x\right)$, whence
$\mathcal{A}(x,o)=\nabla_{x}\mathbb{E}_{Y|x,o}[Y]=\tfrac{\alpha-\beta\gamma}{1+\gamma^{2}}$.

##  Appendix C Simulated study of statins and atherosclerotic cardiovascular
disease

The model is given by

$\displaystyle A$ $\displaystyle=U_{A}$ $\displaystyle L$
$\displaystyle=\theta^{L}_{0}A+\theta^{L}_{1}+e^{\theta^{L}_{2}}U_{L}$
$\displaystyle F$
$\displaystyle=\text{sigmoid}\left(\theta^{F}_{0}+\theta^{F}_{1}(A+\theta^{F}_{2})+\theta^{F}_{3}A^{2}+e^{\theta^{F}_{4}}U_{F}\right)$
$\displaystyle D$
$\displaystyle=\text{sigmoid}\left(\theta^{D}_{0}+\theta^{D}_{1}L+\theta^{D}_{2}A+\theta^{D}_{3}A^{2}+e^{\theta^{D}_{4}}u_{D}\right)$
$\displaystyle R$
$\displaystyle=\text{sigmoid}\left(\theta^{R}_{0}+\theta^{R}_{1}D+\theta^{R}_{2}\log(A)+\theta^{R}_{3}\log^{2}(A)+\theta^{R}_{4}L+\theta^{R}_{5}\log(A)L+\theta^{R}_{6}F\right)$
$\displaystyle X$
$\displaystyle=\text{sigmoid}\Big{(}\theta^{X}_{1}\text{1}_{0.05\leq
R<0.075}+\theta^{X}_{2}\text{1}_{0.075\leq
R<0.2}+\theta^{X}_{3}\text{1}_{R\geq
0.2}+\theta^{X}_{4}+\theta^{X}_{5}\text{1}_{D\geq 0.5}$ $\displaystyle\hskip
48.36958pt+\theta^{X}_{6}L+\theta^{X}_{7}\text{1}_{L>\log(160)}+\theta^{X}_{8}(A+\theta^{X}_{0})+\theta^{X}_{9}(A+\theta^{X}_{0})^{2}+e^{\theta^{X}_{10}}U_{X}\Big{)}$
$\displaystyle M$
$\displaystyle=L+\theta^{M}_{0}+\theta^{M}_{1}X(\theta^{M}_{2}-L)\text{1}_{L<\log(130)}+e^{\theta^{M}_{3}}U_{M}$
$\displaystyle Y$
$\displaystyle=\text{sigmoid}\Big{(}\theta^{Y}_{3}+\theta^{Y}_{4}X+\theta^{Y}_{5}M+\theta^{Y}_{6}\sqrt{A+\theta^{Y}_{0}}+\theta^{Y}_{7}D+\theta^{Y}_{8}e^{1+R}$
$\displaystyle\hskip
48.36958pt+\theta^{Y}_{9}\left(L+\theta^{Y}_{1}\text{sigmoid}\left(10(L+\theta^{Y}_{2})\right)L^{2}\right)+e^{\theta^{Y}_{10}}U_{Y}\Big{)}$
$\displaystyle H$
$\displaystyle=\text{sigmoid}\left(\theta^{H}_{0}+\theta^{H}_{1}X+\theta^{H}_{2}Y+e^{\theta^{H}_{3}}U_{H}\right),$

with

$\displaystyle\theta^{L}$ $\displaystyle=[0.005,\log(100),\log(0.18)]$
$\displaystyle\theta^{F}$ $\displaystyle=[-5.5,0.05,-20,0.001,\log(1.1)]$
$\displaystyle\theta^{D}$ $\displaystyle=[-4.23,0.03,-0.02,0.0009,\log(1.6)]$
$\displaystyle\theta^{R}$ $\displaystyle=[4.3,3.5,-2.07,0.05,4.09,-1.04,0.01]$
$\displaystyle\theta^{X}$
$\displaystyle=[-30,0.273,1.592,2.461,-3.471,1.39,0.112,0.973,-0.046,0.003,\log(1.7)]$
$\displaystyle\theta^{M}$ $\displaystyle=[0.1,-3.5,5,0]$
$\displaystyle\theta^{Y}$
$\displaystyle=[-39,1.4,-\log(110),-6.25,-0.75,-0.1,0.45,1.75,0.29,0.1,\log(0.9)]$
$\displaystyle\theta^{H}$ $\displaystyle=[-1.7,0.8,1.5,\log(0.5)].$

We take $p(U_{A})$ to be a trapezoidal distribution with bottom base from 40
to 75 and top base from 40 to 60. All other distributions of exogenous random
variables are taken to be standard Gaussians, that is
$U_{L},U_{F},U_{D},U_{X},U_{M},U_{Y},U_{H}\sim\mathcal{N}(0,1)$.


# Exact Algorithms for No-Rainbow Coloring and Phylogenetic Decisiveness

Ghazaleh Parvini  College of Information & Computer Sciences, University of
Massachusetts Amherst, Amherst, MA 01003, USA  David Fernández-Baca
Department of Computer Science, Iowa State University, Ames, IA 50011, USA

###### Abstract

The input to the no-rainbow hypergraph coloring problem is a hypergraph $H$
where every hyperedge has $r$ nodes. The question is whether there exists an
$r$-coloring of the nodes of $H$ such that all $r$ colors are used and there
is no rainbow hyperedge — i.e., no hyperedge uses all $r$ colors. The no-
rainbow hypergraph $r$-coloring problem is known to be NP-complete for $r\geq
3$. The special case of $r=4$ is the complement of the phylogenetic
decisiveness problem. Here we present a deterministic algorithm that solves
the no-rainbow $r$-coloring problem in $O^{*}((r-1)^{(r-1)n/r})$ time and a
randomized algorithm that solves the problem in $O^{*}((\frac{r}{2})^{n})$
time.

##  1 Introduction

The hypergraph no-rainbow $r$-coloring problem asks whether there is an
assignment to one of $r$ colors to each of the nodes in a hypergraph such that
all $r$ colors are used, but no edge has $r$ nodes with distinct colors. The
problem is perhaps the most well known in the class of “surjective constraint
satisfaction problems” [2]. It was recently shown to be NP-complete [10], even
for $r=3$.

The case $r=4$ is of interest in building phylogenetic trees from incomplete
data. References [9, 6] formulate the question of determining whether such
data suffices to reconstruct a unique tree as a combinatorial problem, termed
“decisiveness.” The decisiveness problem is equivalent to the complement of
the no-rainbow $4$-coloring problem [5].

In [5] we presented exact algorithms for no-rainbow $r$-coloring for $r=3$ and
$r=4$, that run in time $O^{*}(1.89^{n})$ and $O^{*}(2.81^{n})$ time,
respectively111The $O^{*}$-notation is a variant of $O$-notation that ignores
polynomial factors [3].. Here we present more efficient randomized and
deterministic algorithms based on local search. We note that local search has
been applied earlier to obtain exact algorithms for the $k$-SAT problem [7, 4]
(see also [8]).

##  2 Preliminaries

We write $[r]$ to denote the set $\\{1,2,\dots,r\\}$, for some positive
integer $r$.

A _hypergraph_ $H$ is a pair $H=(X,E)$, where $X$ is a set of elements called
_nodes_ or _vertices_ , and $E$ is a set of non-empty subsets of $X$ called
_hyperedges_ or _edges_ [1]. Two nodes $u,v\in V$ are _neighbors_ if
$\\{u,v\\}\subseteq e$, for some $e\in E$. A hypergraph $H=(X,E)$ is _$r$
-uniform_, for some integer $r>0$, if each hyperedge of $H$ contains exactly
$r$ nodes.

Let $H=(X,E)$ be a hypergraph and $r$ be a positive integer. An _$r$
-coloring_ of $H$ is a mapping $c:X\rightarrow[r]$. For node $v\in X$, $c(v)$
is the _color_ of $v$. Throughout this paper, $r$-colorings are assumed to be
_surjective_ ; that is, for each $i\in[r]$, there is at least one node $v\in
X$ such that $c(v)=i$. The _Hamming distance_ between two colorings $c$ and
$c^{\prime}$ of $H$, denoted $d(c,c^{\prime})$, equals the number of nodes
$v\in X$ such that $c(v)\neq c^{\prime}(v)$. Two colorings $c$ and
$c^{\prime}$ of $H$ _agree_ on set $F\subseteq X$ if $c(v)=c^{\prime}(v)$ for
all $v\in F$.

Let $c$ be an $r$-coloring of $H$. Edge $e\in E$ is a _rainbow edge with
respect to $c$_ if, for each $i\in[r]$, there is at least one $v\in e$ such
that $c(v)=i$. Coloring $c$ _induces a rainbow edge_ if there is an edge $e\in
E$ such that $e$ is a rainbow edge with respect to $c$. A _no-rainbow
$r$-coloring_ of $H$ is a surjective $r$-coloring $c$ of $H$ such that $c$
induces no rainbow edge. Given an $r$-uniform hypergraph $H=(X,E)$, the _no-
rainbow $r$-coloring problem_ asks whether $H$ has a no-rainbow $r$-coloring
[2].

##  3 Deterministic Local Search

Our deterministic algorithm for the no-rainbow $r$-coloring (Algorithm 1)
consists of a procedure DetNRC that conducts multiple searches for a no-
rainbow coloring, starting from different colorings, by invoking a procedure
DetLocalSearch. DetNRC takes as its argument an $r$-uniform hypergraph
$H=(X,E)$ and returns $1$ (true) if $H$ has a no-rainbow $r$-coloring, and 0
(false) otherwise222For simplicity, our algorithm only returns true or false.
The algorithm can be easily modified to return a no-rainbow coloring, if one
exists.. To explain our algorithm in detail, we need to introduce some
concepts.

_DetNRC(_$H$_)_:

foreach  _initial candidate pair $(c,F)$_ do

if  _DetLocalSearch(_$H,(c,F),\frac{(r-1)n}{r}$_)_ then return $1$

return $0$

_DetLocalSearch(_$H,(c,F),g$_)_:

if  _$g=0$ and $c$ induces a rainbow edge in $H$_ then return $0$

if  _$c$ induces a rainbow edge $e$ such that $e\subseteq F$_ then return $0$

if  _$c$ is a no-rainbow coloring_ then return $1$

if  _$|e\cap F|\neq r-1$ for every $e\in E$_ then return $1$

Choose any rainbow hyperedge $e\in E$ such that $|e\cap F|=r-1$

Let $v$ be the node of $e$ such that $v\not\in F$

foreach  _$j\in[r]\setminus\\{c(v)\\}$_ do

Let $c^{\prime}$ be the coloring of $H$ where $c^{\prime}(v)=j$ and
$c^{\prime}(u)=c(u)$ for all $u\in X\setminus\\{v\\}$

if  _$\textnormal{{DetLocalSearch}}(H,(c^{\prime},F\cup\\{v\\}),g-1)$_ then
return $1$

return $0$

Algorithm 1 A deterministic algorithm for no-rainbow $r$-coloring.

A _candidate pair_ for $H=(X,E)$ is a pair $(c,F)$ where $c$ is an
$r$-coloring of $H$ and $F$ is a subset of $X$ such that for each $i\in[r]$
there is a node $v\in F$ such that $c(v)=i$ (thus, $c$ is surjective and
$|F|\geq r$). Once a node is added to set $F$, its color is not allowed to
change. We refer to the nodes in $F$ as _frozen_ nodes.

Starting from an initial candidate pair, DetLocalSearch performs a series of
steps, where each step takes us from the current candidate pair $(c,F)$ to a
new candidate pair $(c^{\prime},F^{\prime})$, where $c^{\prime}$ is obtained
by changing the color of some node $v\not\in F$ to one of the other $r-1$
colors and where $F^{\prime}=F\cup\\{v\\}$. A step may lead to a dead end — a
candidate pair $(c,F)$ where $c$ induces a rainbow edge $e$ and every node in
$e$ is frozen. In this case, we backtrack and try another candidate pair.

To limit the running time of DetLocalSearch, we put a bound on the _search
radius_ to be explored; that is, on the maximum Hamming distance allowed
between the initial coloring and any coloring encountered during the search.
To guarantee correctness, we must choose the starting candidate pairs and the
search radius appropriately.

An _initial candidate pair_ for $H=(X,E)$ is a candidate pair $(c,F)$ where
$|F|=r$ and every node in $X\setminus F$ has the same color (which, obviously,
must be the same as that of one of the nodes in $F$. Note that the number of
initial candidate pairs is ${n\choose r}\cdot r$, which is polynomial for
fixed $r$.

######  Lemma 1.

Let $H=(X,E)$ be an $r$-uniform hypergraph with $n$ nodes. Suppose $H$ has a
no-rainbow $r$-coloring $c^{*}$. Then, there exists an initial candidate pair
$(c,F)$ such that $c$ and $c^{*}$ agree on $F$ and
$d(c,c^{*})\leq\frac{(r-1)n}{r}$.

###### Proof.

Let $(c,F)$ be the initial candidate pair constructed as follows. Set $F$
consists of any $r$ nodes in $X$ that are assigned distinct colors by $c^{*}$.
For each $v\in F$, let $c(v)=c^{*}(v)$. For each $v\in X\setminus F$, let
$c(v)=i$, where color $i\in[r]$ is chosen arbitrarily among the colors used at
least $\frac{n}{r}$ times by $c^{*}$. Then, $(c,F)$ is an initial candidate
pair with $d(c,c^{*})\leq(r-1)n/r$. ∎

The previous lemma puts a limit on the search radius that must be explored,
provided we start with the right candidate pair. In fact, during the search,
we may find candidate pairs $(c,F)$ where $c$ induces a rainbow edge, but yet
we can immediately conclude that a rainbow coloring exists. The next lemma
characterizes this situation.

######  Lemma 2.

Let $H=(X,E)$ be an $r$-uniform hypergraph. Suppose that $(c,F)$ is a
candidate pair for $H$ such that there does not exist a rainbow edge $e$ where
$e\subseteq F$. If $|e\cap F|\neq r-1$ for each edge $e\in E$, then $H$ admits
a no-rainbow $r$-coloring that agrees with $c$ on $F$.

###### Proof.

Let $c^{*}$ be the $r$-coloring of $H$ where $c^{*}(v)=c(v)$ for all $v\in F$
and $c^{*}(v)=1$ for all $v\in X\setminus F$. Then $c^{*}$ is a no-rainbow
$r$-coloring of $H$ because (a) $c^{*}$ is surjective, (b) no edge $e\in E$
such that $|e\cap F|=r$ is a rainbow edge, and (c) for each $e\in E$ such that
$|e\cap F|\leq r-2$, $c$ assigns the same color to at least two nodes in $e$.
∎

Let $(c,F)$ be a candidate pair for a hypergraph $H=(X,E)$. We say that
$(c,F)$ is _within $g$ of a no-rainbow coloring_, for some integer $g\geq 0$,
if there exists a candidate pair $(c^{\prime},F^{\prime})$ where $F\subseteq
F^{\prime}$ such that $c$ and $c^{\prime}$ agree on $F$, $d(c,c^{*})\leq g$,
and either $c^{\prime}$ is a no-rainbow $r$-coloring of $H$ or $|e\cap
F^{\prime}|\neq r-1$ for each edge $e\in E$.

######  Lemma 3.

Let $H$ be an $r$-uniform hypergraph, $(c,F)$ be a candidate pair for $H$ and
$g$ be a positive integer. Then, DetLocalSearch(_$H,(c,F),g$_) returns $1$ if
$(c,F)$ is within $g$ of a no-rainbow coloring, and $0$ otherwise. The running
time of DetLocalSearch is $O^{*}((r-1)^{g})$.

###### Proof.

The first four if statements of DetLocalSearch successively consider the base
cases. If $g=0$ and $c$ induces a rainbow edge, no further color changes are
allowed, so the answer must be $0$. If $c$ induces a rainbow edge $e$ such
that $e\subseteq F$, then we are not allowed to change any of the colors in
$e$, so the answer must again be $0$. If $c$ is a no-rainbow coloring, then we
must obviously return $1$. If $|e\cap F|\neq r-1$ for every edge $e$, then, by
Lemma 2, $H$ has a no-rainbow $r$-coloring, so DetLocalSearch returns $1$.

For the recursion, suppose $(c,F)$ is within $g$ of a no-rainbow coloring and
that $e$ is any rainbow edge such that $|e\cap F|=r-1$. Let $v$ be the single
unfrozen node in $e-F$. Then, there are only $r-1$ possible colors for $v$. By
setting $v$ to one of these alternative colors and adding $v$ to $F$, we
obtain a candidate pair that is within $g-1$ of a no-rainbow coloring.

The running time of DetLocalSearch is $O^{*}((r-1)^{g})$, since the depth of
its recursion tree is $g$, each node has $0$ or $r-1$ children, and the work
per node is polynomial. ∎

######  Theorem 1.

Given an $n$-node $r$-uniform hypergraph $H$, algorithm DetNRC returns $1$ if
$H$ admits a no-rainbow, and $0$ otherwise. The worst-case running time of
DetNRC is $O^{*}((r-1)^{(r-1)n/r})$. In particular, the worst-case running
time is $O^{*}(1.59^{n})$ for $r=3$ and $O^{*}(2.28^{n})$ for $r=4$.

###### Proof.

Suppose $H$ has a no-rainbow $r$-coloring $c^{*}$. Then, by Lemma 1, one of
the initial candidate pairs is within $(r-1)n/r$ of $c^{*}$. Since
DetLocalSearch tries all possible initial candidates, DetNRC is correct.

By Lemma 3, since DetLocalSearch is invoked with search radius
$g=\frac{(r-1)n}{n}$, each call to DetLocalSearch $O^{*}((r-1)^{(r-1)n/r})$
time. Since DetNRC calls DetLocalSearch a polynomial number of times, the
claimed running time follows. ∎

##  4 Randomized Algorithm

Our randomized algorithm for the no-rainbow $r$-coloring problem (Algorithm 2)
consists of a procedure RandNRC whose outer for loop conducts a series of
trials. Each trial invokes a randomized procedure RandLocalSearch to search
for a no-rainbow $r$-coloring starting from different candidate pairs. We will
prove that, given an $n$-node $r$-uniform hypergraph $H$ and a number
$\alpha>1$, RandNRC(_$H,\alpha$_) finds a no-rainbow $r$-coloring of $H$, if
such a coloring exists, with probability at least $1-\frac{1}{e^{\alpha}}$.

_RandNRC(_$H,\alpha$_)_:

for  _$j=1$ to $\alpha\cdot(\frac{r}{2})^{n}$_ do

foreach  _$F\subseteq X$ such that $|F|=r$_ do

Let $c$ be the coloring obtained by assigning distinct colors to the nodes of
$F$ and assigning each node in $X\setminus F$ a color chosen uniformly at
random from $[r]$

if  _RandLocalSearch(_$H,c,F$_)_ then return $1$

return $0$

_RandLocalSearch(_$H,c,F$_)_:

for  _$i=1$ to $n-r$_ do

if  _$c$ induces a rainbow edge such that $e\subseteq F$_ then return $0$

if  _$c$ is a no-rainbow coloring_ then return $1$

if  _$|e\cap F|\neq r-1$ for every edge $e$_ then return $1$

Pick an arbitrary rainbow hyperedge $e$ such that $|e\cap F|=r-1$

Let $v$ be the node of $e$ such that $v\not\in F$

Pick a color $j\in[r]\setminus\\{c(v)\\}$ uniformly at random

$c(v)=j$

$F=F\cup\\{v\\}$

return $0$

Algorithm 2 A randomized algorithm for no-rainbow $r$-coloring.

We begin by bounding the probability that RandLocalSearch finds a no-rainbow
$r$-coloring, provided we start it at the right point.

######  Lemma 4.

Let $H=(X,E)$ be an $n$-node $r$-uniform hypergraph. Suppose $(c,F)$ is a
candidate pair such that $|F|=4$ and for each node $v\in X\setminus F$, $c(v)$
is chosen at random. If $H$ has a no-rainbow $r$-coloring that agrees with $c$
on $F$, then, with probability at least $\left(\frac{2}{r}\right)^{n}$,
RandLocalSearch(_$H,c,F$_) will find a no-rainbow $r$-coloring for $H$.

###### Proof.

Suppose $H$ has a no-rainbow $r$-coloring that agrees with $c$ on $F$.
RandLocalSearch starts from $(c,F)$ and at each step determines that either
(1) it has reached a dead end, because there exists a rainbow edge all whose
nodes are frozen, (2) $c$ is a no-rainbow $r$-coloring, (3) $(c,F)$ satisfies
the conditions of Lemma 2 and, therefore $H$ admits a no-rainbow $r$-coloring
or (4) there is a rainbow edge, say $e$, with exactly one unfrozen node. In
case (1), we say that RandLocalSearch _fails_ , while in cases (2) and (3) it
_succeeds_. In case (4) RandLocalSearch, recolors the unfrozen node, say $v$,
in $e$ with a randomly chosen color, adds $v$ to $F$ and continues the search.

The probability that RandLocalSearch succeeds is the probability that it
reaches cases (2) or (3). Let $(c^{*},F^{*})$ be a candidate pair such that
$c^{*}$ agrees with $c$ on $F$ and $F\subseteq F^{*}$ such that either $c^{*}$
is a rainbow coloring or $(c^{*},F^{*})$ satisfies the conditions of Lemma 2.
By assumption, $(c^{*},F^{*})$ must exist. If at some iteration of
RandLocalSearch we have $(c,F)=(c^{*},F^{*})$, then RandLocalSearch succeeds.

Let $A_{k}$ be the event that the initial assignment $c$ disagrees with
$c^{*}$ on exactly $k$ vertices. By assumption, $c$ agrees with $c^{*}$ on at
least $r$ nodes. Thus, letting $m=n-r$, we have

$\Pr[A_{k}]={m\choose k}(r-1)^{k}r^{-m}.$

The probability $p$ that RandLocalSearch succeeds is at least the probability
that each iteration reduces the distance between $c$ and $c^{*}$ by 1, until
$c=c^{*}$. The probability of getting closer to $c^{*}$ at each step is
$\frac{1}{r-1}$. Thus, we have

$p\geq\sum_{k=0}^{m}{m\choose
k}(r-1)^{k}r^{-m}\left(\frac{1}{r-1}\right)^{k}=r^{-m}\sum_{k=0}^{m}{m\choose
k}=\left(\frac{2}{r}\right)^{m}\geq\left(\frac{2}{r}\right)^{n},$

as claimed. ∎

######  Theorem 2.

Given an $r$-uniform hypergraph $H=(X,E)$ and a number $\alpha>1$, RandNRC
finds a no-rainbow $r$-coloring of $H$, if one exists, with probability at
least $1-\frac{1}{e^{\alpha}}$. The worst-case running time of RandNRC is
$O^{*}((\frac{r}{2})^{n})$. In particular, for $r=3$ and $r=4$, the worst-case
running times are $O^{*}(1.5^{n})$ and $O^{*}(2^{n})$, respectively.

###### Proof.

If $H$ has a no-rainbow $r$-coloring, there must exist a set $F\subseteq X$
with $|F|=r$ such that the nodes in $F$ have $r$ different colors. Thus, by
Lemma 4, at least one of the calls to RandLocalSearch has probability at least
$\left(\frac{2}{r}\right)^{n}$ of success. It follows from standard results
(see, e.g., [8, p. 151]) that by conducting $\alpha\cdot(\frac{r}{2})^{n}$
trials of RandNRC’s outer for loop, the probability of success is at least
$1-\frac{1}{e^{\alpha}}$. The time bound follows by noting that the running
time of RandLocalSearch is polynomial and, for each $j$, the number of
executions of RandNRC’s inner for loop is polynomial for fixed $r$. ∎

## References

  * [1] Claude Berge.  Hypergraphs: combinatorics of finite sets, volume 45 of North-Holland Mathematical Library.  Elsevier, 1989.
  * [2] Manuel Bodirsky, Jan Kára, and Barnaby Martin.  The complexity of surjective homomorphism problems—a survey.  Discrete Applied Mathematics, 160(12):1680–1690, 2012.
  * [3] Fedor V. Fomin and Dieter Kratsch.  Exact Exponential Algorithms.  Springer, Berlin, Heidelberg, 2010.
  * [4] Robin A Moser and Dominik Scheder.  A full derandomization of Schöning’s $k$-SAT algorithm.  In Proceedings of the 43rd Annual ACM Symposium on Theory of Computing, pages 245–252, 2011.
  * [5] Ghazaleh Parvini, Katherine Braught, and David Fernández-Baca.  Checking phylogenetic decisiveness in theory and in practice.  In Zhipeng Cai, Ion I. Mandoiu, Giri Narasimhan, Pavel Skums, and Xuan Guo, editors, Bioinformatics Research and Applications - 16th International Symposium, ISBRA 2020, Moscow, Russia, December 1-4, 2020, Proceedings, volume 12304 of Lecture Notes in Computer Science, pages 189–202. Springer, 2020.
  * [6] Michael J. Sanderson, Michelle M. McMahon, and Mike Steel.  Terraces in phylogenetic tree space.  Science, 333(6041):448–450, 2011.
  * [7] Uwe Schöning.  A probabilistic algorithm for k-sat and constraint satisfaction problems.  In 40th Annual Symposium on Foundations of Computer Science (Cat. No. 99CB37039), pages 410–414. IEEE, 1999.
  * [8] Uwe Schöning and Jacobo Torán.  The Satisfiability Problem: Algorithms and Analyses, volume 3 of Mathematik für Anwendungen.  Lehmanns Media, Berlin, 2013.
  * [9] Mike Steel and Michael J. Sanderson.  Characterizing phylogenetically decisive taxon coverage.  Applied Mathematics Letters, 23(1):82–86, 2010.
  * [10] Dmitriy Zhuk.  No-rainbow problem is NP-hard, March 2020.  arXiv preprint https://arxiv.org/abs/2003.11764.