# Stochastic Causal Programming for Bounding Treatment Effects

Kirtan Padh  Helmholtz AI, Munich  Technical University of Munich  Jakob
Zeitler  University College London  David Watson  University College London
Matt Kusner  University College London  Ricardo Silva  University College
London  Niki Kilbertus  Helmholtz AI, Munich  Technical University of Munich

###### Abstract

Causal effect estimation is important for numerous tasks in the natural and
social sciences. However, identifying effects is impossible from observational
data without making strong, often untestable assumptions. We consider
algorithms for the _partial_ identification problem, bounding the effects of
multivariate, continuous treatments over multiple possible causal models when
unmeasured confounding makes identification impossible. We propose a framework
where observable evidence is matched to the implications of constraints
encoded in a causal model by norm-based criteria. This generalizes classical
approaches based on generative models. Casting causal effects as objective
functions in a constrained optimization problem, we combine flexible learning
algorithms with Monte Carlo methods to implement a family of solutions under
the name of _stochastic causal programming_. In particular, we present ways by
which such constrained optimization problems can be parameterized without
likelihood functions for the causal or the observed data model, reducing the
computational and statistical complexity of the task.

##  1 Introduction

Estimating causal effects is a key goal of scientific inquiry, enabling better
decision making in medicine (Bica et al., , 2021), economics (Chernozhukov et
al., , 2018), epidemiology (Vandenbroucke et al., , 2016), and beyond (Pearl,
, 2009; Imbens and Rubin, , 2015). The gold standard for estimating these
effects is the randomized controlled trial (RCT): randomly assign individuals
to treatment and control groups (e.g., vaccine or placebo) and observe
outcomes (e.g., antibody levels). These experiments work because randomization
removes _confounding_ , which occurs when pre-treatment covariates drive both
treatment propensities and potential outcomes. However, in many cases it is
physically, logistically, or ethically impossible to conduct RCTs. In such
cases, causal effects can only be identified under structural assumptions
(Pearl, , 2009; Hernán and Robins, , 2020).

When unobservable confounders are present, techniques such as instrumental
variable (IV) models (Hernán and Robins, , 2020, Ch. 16), the front-door
criterion (Pearl, , 2009, Ch. 3) and proxies (Tchetgen et al., , 2022) may be
applicable. Even then, identification requires further assumptions, not just
about the graphical structure but about the functional form of causal
associations, e.g., monotonicity (Angrist and Imbens, , 1995), additivity
(Hartford et al., , 2017; Singh et al., , 2019; Muandet et al., , 2020;
Bennett et al., , 2019) or conditions such as completeness (Tchetgen et al., ,
2022), which are not immediately intuitive as they refer to unspecified hidden
variables. When such assumptions fail, effects may still be _partially
identifiable_ —i.e., bounded with respect to the set of models consistent with
the data (Manski, , 1990). This task has not received nearly as much attention
as the point estimation problem, despite promising early work in this area
(Chickering and Pearl, , 1996; Balke and Pearl, , 1997). Lately the topic has
sparked some interest in reinforcement learning (Kallus and Zhou, , 2020;
Zhang and Bareinboim, , 2020), algorithmic fairness (Wu et al., 2019a, ; Wu et
al., 2019b, ), and other settings (Gunsilius, , 2019; Zhang and Bareinboim, ,
2021; Hu et al., , 2021). Xia et al.,  (2021) outline a procedure for bounding
causal effects with neural networks but restrict their focus to discrete, low-
dimensional data. Duarte et al.,  (2021); Zhang et al.,  (2021) reduce the
bounding problem to a polynomial program, but both assume that variables are
discrete and Zhang et al.,  (2021) further assumes that variables are finite.
In a recent paper, Kilbertus et al.,  (2020) propose a method for computing
causal bounds in IV models with continuous treatments using gradient descent
and Monte Carlo integration. Their procedure is limited to univariate settings
(see Appendix G for further limitations).

We extend this work to higher dimensions and different causal structures,
replacing the copula formulation of Kilbertus et al.,  (2020) with a more
generic parameterization that decouples model fitting from changes due to the
unobserved confounder $U$.

This parameterization is related to recent literature on (conditional)
normalizing flows (Papamakarios et al., , 2021), where base distributions are
converted into targets through some diffeomorphic transformation(s). Whereas
the goal in flow models is typically generative, the sampling step in our
method is a means toward the end of causal effect bounding.

We make the following contributions: (1) We propose a generic, modular form of
the effect bounding problem compatible with a wide range of graphical
structures, function classes, optimization procedures, and distance measures,
as well as a systematic way of simplifying the computation. (2) We derive an
efficient solution to the problem via constraint sampling and automatic
differentiation, extending the work of Kilbertus et al.,  (2020) to data of
arbitrary dimension. (3) We illustrate our method on synthetic and semi-
synthetic datasets, where simulations confirm that our approach computes valid
bounds even in complex settings where common assumptions fail. We first
present our ideas in a general framework in Sections 2, along with
illustrative examples. In Section 3, we introduce a concrete instantiation of
the framework, with explicit assumptions about the model space and
optimization procedure. Following a number of experiments in Section 4,
Section 5 concludes.

##  2 Setup

We first formulate the general setup and then provide two concrete running
examples based on the IV model and a “leaky mediation” setting. Assume we
observe data from causal model $S^{\star}$, with density function
$p_{S^{\star}}$, including a (possibly multivariate) continuous treatment
$X\in\mathbb{R}^{p}$ and outcome $Y\in\mathbb{R}$ among the observed
variables. Using Pearl’s _do_ notation Pearl,  (2009), the estimand of
interest is
$o_{x^{\star}}(S^{\star}):=\operatorname{\mathbb{E}}_{S^{\star}}[Y\,|\,do(X=x^{\star})]$,
typically not identified. Our goal will be to minimize/maximize
$o_{x^{\star}}(S^{\star})$ over all _admissible_ causal models $S$ among a
(uncountable) model class $\mathcal{S}$. We characterize admissibility via two
types of constraints. First, the missing edges of the assumed directed acyclic
graph (DAG) encode conditional independencies between variables. Further
assumptions on the functional form of structural equations are required for
non-vacuous solutions (Kilbertus et al., , 2020; Gunsilius, , 2021). We refer
to the constraints implied by graphical and functional assumptions jointly as
_structural constraints_. Second, we have observational data from $S^{\star}$.
Denoting an estimate of the data distribution $p_{S^{\star}}$ by $\hat{p}$,
the distribution $p_{S}$ entailed by an admissible causal model $S$ should
replicate this distribution. We measure this _data constraint_ via a distance
function $\operatorname{\mathrm{dist}}(p_{S},\hat{p})$ into $\mathbb{R}^{+}$
that measures the (estimated) discrepancy between model $S\in\mathcal{S}$ and
the ground truth model $S^{\star}$ via observations.

Hence, we arrive at the following general problem setting for computing the
minimum/maximum causal effect among all models $\mathcal{S}$ that are
$(\operatorname{\mathrm{dist}},\epsilon)$-compatible with the observed data:

$\displaystyle\underset{S\in\mathcal{S}}{\text{min / max}}\quad$
$\displaystyle o_{x^{\star}}(S)$ [obj] (1) subject to
$\displaystyle\operatorname{\mathrm{dist}}(p_{S},\hat{p})\leq\epsilon\>,$
[c-data] (2) structural constraints [c-struct] (3)

Note that, unlike most optimization objectives in machine learning, the
objective in Equation 1 depends on the data _only_ via (functionals of) a
fitted distribution $\hat{p}$. This suggests an interpretation of our task as
decoupling model selection from model fitting: the former is a search over an
_equivalence class_ of causal models that (approximately) match the observable
distribution. The latter is in the realm of standard machine learning, which
produces (functionals of) $\hat{p}$ regardless of causal assumptions. We
believe this to be a fruitful viewpoint for a large variety of problems in
causal inference that require us to simultaneously accommodate different
qualitative structural assumptions and provide a good fit to observed data.
Instead of conflating regularization and causal specification when no unique
causal model explains the data (implying likelihood functions full of high-
dimensional plateaus), our formulation considers all models in reasonable
agreement with the estimated observable distribution given by any black-box
algorithm. Importantly, as we show later, we can pick distance functions such
that no full density estimation problem has to be solved.

Consider the instrumental variable model in Figure 1a. Following the
_structural causal model_ (SCM) notation of Pearl (Pearl, , 2009), we observe
instruments $Z\in\mathcal{Z}\subset\mathbb{R}^{q}$, treatments
$X\in\mathcal{X}\subset\mathbb{R}^{p}$, and outcomes
$Y\in\mathcal{Y}\subset\mathbb{R}$ with a potentially high-dimensional
confounded pair of _background_ variables $U_{X}$ and $U_{Y}$. Here we assume
(A1) $X\operatorname{\not\\!\perp\\!\\!\\!\perp}Z$, (A2)
$Z\operatorname{\perp\\!\\!\\!\perp}\\{U_{X},U_{Y}\\}$, and (A3)
$Z\operatorname{\perp\\!\\!\\!\perp}Y\,|\,\\{X,U_{X},U_{Y}\\}$. We denote the
i.i.d. observations that make up $\hat{p}$ by
$\mathcal{D}=\\{x_{i},y_{i},z_{i}\\}_{i=1}^{n}$. Instead of the common
additive noise or linearity assumptions, we consider a larger class of
potentially _non-linear_ , _non-additive_ functions
$f:\mathcal{X}\times\mathcal{U}_{Y}\to\mathcal{Y},g:\mathcal{Z}\times\mathcal{U}_{X}\to\mathcal{X}$
that we will later encode in $\mathcal{S}$. We propose to avoid full density
estimation for [c-data] and instead only estimate $X\,|\,Z$ and the first two
moments of $Y\,|\,\\{X,Z\\}$.

In the SCM of the leaky mediation setting in Figure 1d., a second confounder
between $M$ and $Y$ is added to the standard non-parametrically identified
mediation example for front-door adjustment (Pearl, , 2009, Ch. 3), rendering
the causal effect of $X$ on $Y$ unidentifiable. Here, we have the structural
constraints (B1) $X\operatorname{\perp\\!\\!\\!\perp}Y\,|\,\\{M,U_{Y}\\}$ and
(B2) $U_{X}\operatorname{\perp\\!\\!\\!\perp}U_{M}$, and again allow for _non-
linear_ , _non-additive_ functions
$f:\mathcal{M}\times\mathcal{U}_{Y}\to\mathcal{Y},g:\mathcal{X}\times\mathcal{U}_{M}\to\mathcal{M}$.

Requirements. The key difficulty in operationalizing the optimization in
Equation 1 is the parameterization of $\mathcal{S}$. Our design choices are
guided by efficiency of  (i) the (constrained) optimization of
$o_{x^{\star}}(S)$,  (ii) the evaluation of
$\operatorname{\mathrm{dist}}(p_{S},\hat{p})$,  (iii) restricting
$\mathcal{S}$ to causal models satisfying [c-struct].  Moreover, we explicitly
aim for methods that scale to large sample sizes $n$ as well as multivariate
treatments $(p\\!>\\!1)$ and instruments/mediators $(q\\!>\\!1)$.

##  3 Method

We start this section by first describing the suggested function space
representation (Section 3.1) of SCMs, and how to avoid unnecessary causal
parameters by exploiting some structural signatures (Section 3.2). This is
followed by a detailed algorithmic description of the optimization procedure.

###  3.1 Response function framework

Specifying a causal model $S$, such as the one depicted in Figure 1a.,
involves functions such as $f$ and $g$, as well as the distribution of the
potentially infinite dimensional confounders $U=\\{U_{X},U_{Y}\\}$. To tackle
this complexity, we make use of the so-called _response function_ framework
(Balke and Pearl, , 1994; Kilbertus et al., , 2020). We call
$f_{u}(\cdot):=f(\cdot,u):\mathcal{X}\to\mathcal{Y}$ the response function for
the fixed value $U=u$ of the background variable $U$. Note that in this case,
$f_{u}$ indeed encodes the direct causal effect of $X$ on $Y$. A distribution
over values of $U$ entails a distribution over response functions $f_{u}$.
Hence, instead of explicitly modeling $U$ and $f$ separately as part of a
causal model $S$, we can directly encode them jointly via a distribution over
response functions.

Without restrictions on the functional dependence of $X$ and $Y$ on $U$, one
cannot obtain non-trivial bounds (Gunsilius, , 2021, 2019). As in general we
have little information about the dimensionality and distribution of $U$, we
argue that it is more practical to work with the equivalent assumptions on the
function space and distribution of response functions $f_{u}$. Specifically,
we choose a family of response functions
$\mathcal{F}\subset\\{f:\mathcal{X}\to\mathcal{Y}\\}$ that captures plausible
direct causal effects from $X$ to $Y$ under all possible values of the
confounding variable(s). We then assume a parameterized family of
distributions over $\mathcal{F}$, denoted by
$\\{p_{\eta}^{\mathcal{F}}\,|\,\eta\in\mathbb{R}^{d}\\}$. For concreteness, we
propose parameterized response functions as linear combinations of a set of
non-linear basis functions
$\\{\psi_{k}:\mathcal{X}\to\mathcal{Y}\\}_{k=1}^{K}$ for some
$K\in\mathbb{N}$:

$\mathcal{F}:=\Bigl{\\{}f_{\theta}:=\sum_{k=1}^{K}\theta_{k}\psi_{k}\,\Big{|}\,\theta_{k}\in\mathbb{R}\Bigr{\\}}\>.$
(4)

Hence, once $\mathcal{F}$ is fixed, $p_{\eta}^{\mathcal{F}}$ is fully
described as a distribution over $\theta\in\mathbb{R}^{K}$. For ease of
notation, we often drop the dependence on $\mathcal{F}$ and summarize
$\psi:=(\psi_{1},\ldots,\psi_{K})$. We could fit $\psi_{k}$ adaptively with
almost no modification to the algorithm presented in the sequel (as long as
they are differentiable), but we will choose to keep the basis functions fixed
to demonstrate in the experiments the trade-off between other sources of
causal knowledge and the informativeness of the bounds.

A naïve interpretation of Equation (4), combined with a fixed choice of basis
functions that we adopt for the sake of illustration, is that it may introduce
unwanted constraints in the causal model and may lead to “invalid” bounds. We
address these possible misunderstandings in Appendix A.

###  3.2 Graphical reduction

We reduce a given graphical structure to simplify computation and remove
unneeded assumptions. Figure 1 illustrates the following reduction steps for
IV and leaky mediator examples.

  1. a.

Imagine we are given a graph that explicitly includes the latent background
variables $U$ of a _structural causal model_ (Pearl, , 2009), where each
observed variable is a deterministic function of its direct causes. We express
unmeasured confounding with bidirected edges among the background variables.
This forms a model of _marginal independence_ (Richardson, , 2003) among $U$
variables, which can be hard to parameterize.

  2. b.

To simplify computation, we convert as many bidirected edges to directed edges
as possible while preserving the same structural constraints (Drton and
Richardson, , 2008). This allows for a simplified parameterization. For
instance, for the leaky mediator we have
$p(u_{x},u_{m},u_{y})=p(u_{x})\,p(u_{m})\,p(u_{y}\,|\,u_{x},u_{m})$.

  3. c.

We can remove unnecessary assumptions about structural equations when
$p(v\,|\,do(\mathrm{pa}(v))$, the distribution of some variable $V$ under
intervention on its parents $\mathrm{pa}(V)$, is _identifiable_. Instead of
treating $U_{V}$ as structural background variables with an unknown
distribution and dimensionality, we convert them into “non-causal” random
variables $N_{V}$ of the same dimension as $V$, with a fixed distribution so
that $V$ is a deterministic function of its parents, and which is invertible
on $N_{V}$. For instance, we can model $p(x\,|\,z,u_{x})$ in the IV model as
$X=F_{X\,|\,Z}^{-1}(N_{X})$, where $X$ is one-dimensional and
$F_{X\,|\,Z}(\cdot)$ is the cdf of $X$ given $Z$ and $N_{X}$ is uniform on
$[0,1]$.

  4. d.

Finally, any remaining background variable will require further causal
assumptions. For that, we rely on recasting the conditional model for each
variable $V$ with a parent $U_{V}$ as being generated by a mixture of
deterministic functions, as explained in the previous section.

All of the above provides an alternative to blackbox models such as those of
Hu et al.,  (2021), which ignore that many of the assumptions about structural
equations can be avoided, and that many factors of the data distribution can
be directly encoded without wasteful Monte Carlo simulations.

0pt  IV setting a.$Z$$X$$Y$$U_{X}$$U_{Y}$ b.$Z$$X$$Y$$U_{X}$$U_{Y}$
c.$Z$$X$$Y$$N_{X}$$U_{Y}$ d.$Z$$X$$Y$$N_{X}$$\eta_{Y}$ 0pt Leaky mediation
$X$$\in\mathbb{R}^{p}$ treatments$M$$\in\mathbb{R}^{q}$
mediators$Y$$\in\mathbb{R}$ outcome$U_{X}$$U_{M}$$U_{Y}$
$X$$M$$Y$$U_{X}$$U_{M}$$U_{Y}$ $X$$M$$Y$$N_{X}$$N_{M}$$U_{Y}$
$X$$M$$Y$$N_{X}$$N_{M}$$\eta_{Y}$

(a) * (b) * Figure 1: We show the general graphical reduction procedure for
modeling latent confounding described in the text for the IV and leaky
mediation setting.

###  3.3 Parameterizing $\mathcal{S}$ and satisfying [c-struct]

Next, we describe how to leverage the response function framework for a useful
parameterization of causal models $\mathcal{S}$. Again, we focus on the IV
case for concreteness; see Appendix B for our analysis of the leaky mediator
setting. Despite recent advances in adding exact (conditional) independencies
[c-struct] as explicit constraints in a smooth gradient-based optimization
(Zheng et al., , 2018), this approach is still unwieldy in our case with
unobserved confounding. Therefore we encode [c-struct] directly in our choice
of model class $\mathcal{S}$ and write for clarity

$\overline{\mathcal{S}}:=\\{p_{\eta}\,|\,\eta\in\mathbb{R}^{d};\text{(A1)-(A3)
satisfied for }(X,Y,Z)\sim p_{\eta}\\}\>.$

Using the parameterization from Equation 4, we may now reformulate (A2) and
(A3) as $Z\operatorname{\perp\\!\\!\\!\perp}\theta$ and
$Z\operatorname{\not\\!\perp\\!\\!\\!\perp}\theta\,|\,X$. Consequently, we
propose the DAG in Figure 1b. as a model for $\theta$ that explicitly encodes
(A2) and (A3) via $d$-separation (Pearl, , 2009). Hence, when basing all
$p_{\eta}\in\overline{\mathcal{S}}$ on this graphical structure, we have
absorbed [c-struct] into the definition of $\overline{\mathcal{S}}$. We think
of $N$ as a source of randomness111For emphasis, we repeat that $N$ is _not_
to be interpreted as causal background variables for $X$: it is merely a
device for parameterizing $p(x~{}|~{}z)$ and $p(\theta~{}|~{}x,z)$, and it is
not meant to have a causal interpretation. that is shared in modeling
$X\,|\,\\{Z,N\\}$ and $\theta\,|\,N$, thereby still allowing for dependencies
between $X$ and $\theta$. Observe that a model for $X\,|\,\\{Z,N\\}$ can be
fit once upfront from observed data and remains fixed throughout the
optimization. Later, we will need to deterministically recover $N$ from $Z$
and $X$, leaving us with any invertible conditional model $h_{Z}(N)$ for which
$h_{Z=z}^{-1}(x)$ exists uniquely for all $z\in\mathcal{Z}$ and
$x\in\mathcal{X}$. For example, conditional normalizing flows composed of
invertible and differentiable (i.e., diffeomorphic) transformations
(Papamakarios et al., , 2021) are flexible candidates that we use in our
experiments. Note that this parameterization implies $d_{N}=p$. We provide all
details of the specific implementation in Appendix E.

Once $h_{Z}$ is fixed, we model $\theta$ with a _partially specified_
multivariate distribution as follows,

$\theta\,|\,N\sim
p_{\eta.\theta}\bigl{(}\cdot;\mu_{\eta_{0}}(N),\Sigma_{\eta_{1}}(N)\bigr{)},$
(5)

with mean and covariance functions $\mu_{\eta_{0}}(N)$ and
$\Sigma_{\eta_{1}}(N)$, both parameterized by our optimization parameters
$\eta=(\eta_{0},\eta_{1})$. Other than that, we make _no further assumptions_
about the shape of $p_{\eta.\theta}$, as they are unnecessary for the
inference problems we set up to solve. To ensure $\Sigma_{\eta_{1}}(N)$ is a
valid covariance matrix, we represent it in terms of its Cholesky factor $L$
and add a small constant $\Omega$ to the diagonal
$\Sigma_{\eta_{1}}=L^{\top}L+\Omega\bm{1}$. Again we can use flexible function
approximators for $\mu_{\eta_{0}}$ and $\Sigma_{\eta_{1}}$ and resort to
neural networks with parameters $\eta_{0}$ and $\eta_{1}$ in our experiments.
The combined parameters $\eta=(\eta_{0},\eta_{1})$ now encode our family of
causal models
$\overline{\mathcal{S}}=\\{p_{\eta}\,|\,\eta\in\mathbb{R}^{d}\\}$, where $Y$
is implicitly given by the random variable $f_{\theta}(X)$ with $X=h_{Z}(N)$,
$N\sim\mathcal{N}(\bm{0},\bm{1})$, and the required moments of $\theta\,|\,N$
as in Equation 5. We remark that the choice of which functionals/moments of
$\mathcal{F}_{\theta}$ to parameterize depends on the information from the
data that we want to use, which in the sequel will be the first and second
moments of $Y~{}|~{}\\{X,Z\\}$.

We now rephrase our main optimization problem as:

$\displaystyle\underset{\eta\in\mathbb{R}^{d}}{\text{min / max}}\quad$
$\displaystyle
o_{x^{\star}}(\eta)=\psi(x^{\star})^{\top}\operatorname{\mathbb{E}}_{N}[\mu_{\eta_{0}}(N)]$
[obj] subject to
$\displaystyle\operatorname{\mathrm{dist}}(p_{\eta},\hat{p})\leq\epsilon\>.$
[c-data]

For the objective equation, we used the law of total expectation. Since $N$
follows a standard Gaussian, we can easily estimate the remaining expectation
from a finite sample. So far, the data have only entered once when we chose
$\overline{\mathcal{S}}$ such that $\hat{p}(X\,|\,Z)$ is matched by all
$p_{\eta}\in\overline{\mathcal{S}}$. We still must enforce the remaining data
matching constraints [c-data].

Algorithm 1 Computing upper or lower bounds on
$\operatorname{\mathbb{E}}[Y\,|\,do(X=x^{\star})]$ in the IV setting.

1:dataset $\mathcal{D}=\\{(z_{i},x_{i},y_{i})\\}_{i=1}^{n}$; constraint
functions $\\{\phi_{l}:\mathcal{Y}\to\mathbb{R}\\}_{l=1}^{L}$; basis functions
$\\{\psi_{k}:\mathcal{X}\to\mathcal{Y}\\}_{k=1}^{K}$; norm $\|\cdot\|$ for
$\operatorname{\mathrm{dist}}$; batchsize for Monte Carlo $B$; number of
support points $M$; tolerance $\epsilon>0$

2:Setup: One-time computations shared for different $x^{\star}$ values

3:Fit (invertible) conditional normalizing flow $X=h_{Z}(N)$ from data
$\mathcal{D}$ for $N\sim\mathcal{N}(0,\bm{1}_{p})$

4:Fit MLPs $\hat{\phi}_{1}:x_{i},z_{i}\to y_{i}$ and
$\hat{\phi}_{2}:x_{i},z_{i}\to y_{i}^{2}$ by minimizing the squared loss from
data $\mathcal{D}$

5:subsample $M$ indices from $[n]$ (uniform, no replacement) $\triangleright$
“support points”, w.l.o.g. use $[M]$

6:Optimization: performed separately for lower and upper bound for each
$x^{\star}$

7:minimize Objective($\eta$) subject to Constraint($\eta$)${}\leq\epsilon$
(see Section E.5)

8:

9:function Objective($\eta$)

10:
$o_{x^{\star}}(\eta)\leftarrow\psi(x^{\star})^{\top}\frac{1}{B}\sum_{j=1}^{B}\mu_{\eta_{0}}(n_{j})$
with $n_{j}\sim\mathcal{N}(0,\bm{1}_{p})$ $\triangleright$ differentiable
w.r.t. $\eta$

11: return $\pm o_{x^{\star}}(\eta)$ $\triangleright$ objective, $\pm$ for
lower/upper bound

12:function Constraint($\eta$)

13: $n_{j}\leftarrow h_{z_{j}}^{-1}(x_{j})$ for $j\in[M]$ $\triangleright$
invert $X\,|\,Z$ model to infer “noises”

14: $A_{1,j}(\eta)\leftarrow\psi(x_{j})^{\top}\mu_{\eta_{0}}(n_{j})$ for
$j\in[M]$ $\triangleright$ moments implied by model

15:
$A_{2,j}(\eta)\leftarrow\psi(x_{j})^{\top}\Bigl{(}\Sigma_{\eta_{1}}(n_{j})+\mu_{\eta_{0}}(n_{j})\mu_{\eta_{0}}(n_{j})^{\top}\Bigr{)}\psi(x_{j})$
for $j\in[M]$

16: $\nu_{l,j}\leftarrow\hat{\phi}_{l}(x_{j},z_{j})-A_{l,j}(\eta)$ for
$l\in\\{1,2\\},j\in[M]$ $\triangleright$ constraint matrix

17: return $\|\nu\|$

###  3.4 Satisfying [c-data]: matching the data

We first factorize $p(Z,X,Y)=p(Y\,|\,X,Z)\,p(X\,|\,Z)\,p(Z)$. A key advantage
of our construction of $\overline{\mathcal{S}}$ is that we only need to match
the factor $p(Y\,|\,X,Z)$ to satisfy [c-data] instead of having to perform
density estimation to match high-dimensional distributions via, e.g.,
f-divergences, integral probability metrics, Stein discrepancy, adversarial
optimization, or variational inference. For us, any efficiently computable
similarity measure $\operatorname{\mathrm{dist}}$ on $\hat{p}(Y\,|\,X,Z)$ and
$f_{\theta}(X)\,|\,\\{X,Z\\}$ can be used to keep candidate models $p_{\eta}$
close to matching the _entire_ observed data distribution. This matching can
be conveniently formulated as matching expectations of the two distributions
when transformed by a fixed set of dictionary functions
$\\{\phi_{l}:\mathcal{Y}\to\mathbb{R}\\}_{l=1}^{L}$. That is, we want the
absolute values of

$\nu_{l}(x,z)=\operatorname{\mathbb{E}}_{y\sim\hat{p}(y\,|\,x,z)}[\phi_{l}(y)]-\operatorname{\mathbb{E}}_{\theta\sim
p_{\eta}(\theta\,|\,x,z)}[\phi_{l}(f_{\theta}(x))]\>$ (6)

to be small for all $l$ and all $(x,z)\in\mathcal{X}\times\mathcal{Z}$. Since
$X,Z$ are continuous (and potentially high-dimensional), this would amount to
an uncountably infinite set of constraints. In high dimensions, a naïve
discretization of $X,Z$ also fails (Kilbertus et al., , 2020). Instead, we aim
to make $|\nu_{l}(x,z)|$ small at a representative, finite set of points in
$\mathcal{X}\times\mathcal{Z}$. Arguably, the most representative set is a
uniformly random subsample of the observed data $\mathcal{D}$ of size $M$. For
notational simplicity, we assume w.l.o.g. that these “support points” are the
first $M$ indices. In other words, defining (by some overloading of notation)
the matrix $\nu\in\mathbb{R}^{L\times M}$ via $\nu_{l,j}=\nu_{l}(x_{j},z_{j})$
(using $j=1,\ldots,M$ to index the selected support points), we obtain a
natural family of overall distances
$\operatorname{\mathrm{dist}}(p_{\eta},\hat{p})=\|\nu\|$ for any (semi-)norm
$\|\cdot\|$ on $\mathbb{R}^{L\times M}$.222We remark that, formally, if we
allow $L,M\to\infty$ the functions $\phi_{l}$ can be chosen such that exact
matching between $\hat{p}$ and $p_{\eta}$ can be approximated arbitrarily
well, e.g., by approximating the characteristic function using $l$ and
matching it everywhere on $\mathcal{X}\times\mathcal{Y}$ using the entry-wise
sup-norm. In our experiments, we specifically consider the entry-wise sup-norm
$\|\cdot\|_{\infty,\infty}$ and the entry-wise 2-norm $\|\cdot\|_{2,2}$. For
$\|\cdot\|_{\infty,\infty}$ we require the absolute value of each entry to be
small, i.e., practically have to enforce $M\cdot L$ constraints, whereas for
$\|\cdot\|_{2,2}$ we get away with a single overall constraint.

The first expectation in Equation 6 is estimable from the observed data via
$\hat{\phi}_{l}(x,z)\approx\operatorname{\mathbb{E}}[\phi_{l}(Y)\,|\,x,z]$,
where $\hat{\phi}_{l}:\mathcal{X}\times\mathcal{Z}\to\mathbb{R}$ are
regression functions mapping
$x,z\to\operatorname{\mathbb{E}}[\phi_{l}(Y)\,|\,x,z]$ that were learned in a
supervised fashion on the training data $\mathcal{D}$. Because $h_{Z}$ was
assumed to be invertible, our modeling choice for $\overline{\mathcal{S}}$
plays nicely with our constraint formulation, in that we can explicitly write

$\theta\,|\,x,z\sim\mathcal{F}_{\theta}\Bigl{(}\theta;\mu_{\eta_{0}}\bigl{(}h_{z}^{-1}(x)\bigr{)},\Sigma_{\eta_{1}}\bigl{(}h_{z}^{-1}(x)\bigr{)}\Bigr{)}\>.$
(7)

For general dictionary functions $\phi_{l}$, we can again estimate the second
expectation in Equation 6 using finite samples from a fully-specified version
of Equation 7. To simplify optimization and modeling, in practice it often
suffices to constrain the first few moments of $Y\,|\,\\{X,Z\\}$, e.g., to set
$\phi_{1}(Y)=Y$ and $\phi_{2}(Y)=Y^{2}$. This task is much more data efficient
than modeling the complete joint distribution $p(X,Y,Z)$, as required, for
example, by generative methods such as GANs (Hu et al., , 2021). This further
simplifies the second expectation in Equation 6 to the following:

$A_{1,j}(\eta)=\psi(x_{j})^{\top}\mu_{\eta_{0}}(n_{j})\>,\;\quad
A_{2,j}(\eta)=\psi(x_{j})^{\top}\bigl{(}\Sigma_{\eta_{1}}(n_{j})+\mu_{\eta_{0}}(n_{j})\mu_{\eta_{0}}(n_{j})^{\top}\bigr{)}\psi(x_{j})\>,$
(8)

where $n_{j}:=h_{z_{j}}^{-1}(x_{j})$ for $j\in[M]$.

Figure 2: Top Left (IV-lin-2d-weak; comparing norms) Our framework easily
allows for different data-matching criteria ($\operatorname{\mathrm{dist}}$).
The choices shown here ($\|\cdot\|_{\infty,\infty}$ and $\|\cdot\|_{2,2}$)
yield comparable and consistent bounds. Top Right (LM-lin-2d-strong; leaky
mediator setting) We get reliable bounds also in the leaky mediator setting
compared to the GAN framework. Bottom Left (IV-lin-2d-strong; comparing
response functions) As expected, neural response functions give wider bounds
than a linear polynomial basis because of being less expressive. Bottom Right
(IV-lin-2d-strong; GAN comparison) Our method reliably yields valid bounds
even under strong confounding, while the GAN framework becomes unstable,
potentially due to difficulties with the bilevel optimization problem.

###  3.5 Solving the optimization

Taking all the steps from previous sections together, we have
$\nu_{l,j}=\hat{\phi}_{l}(x_{j},z_{j})-A_{l,j}(\eta)$ and ultimately arrive at
the following non-convex optimization problem with non-convex constraints

$\displaystyle\underset{\eta\in\mathbb{R}^{d}}{\text{min / max}}\;$
$\displaystyle
o_{x^{\star}}(\eta)=\psi(x^{\star})^{\top}\operatorname{\mathbb{E}}_{N}[\mu_{\eta_{0}}(N)]$
[obj] (9) subject to $\displaystyle\|\nu\|\leq\epsilon$ [c-data]

for which the augmented Lagrangian method with inequality constraints is a
natural choice (Nocedal and Wright, , 2006, Sect. 17.4). We provide a high-
level description in Algorithm 1 and provide specifics in Section E.5. Note
that for the entry-wise sup-norm $\|\cdot\|_{\infty,\infty}$ we aim at
enforcing $M\cdot L$ inequality constraints, one for each entry. While we have
already argued that $L=2$ is a sensible choice, we get to choose a fixed
number of “support datapoints” $M$ regardless of the dimensionality of $X$ to
trade off capturing the observed distribution and stability of the augmented
Lagrangian optimization. For very high-dimensional data, one may take
inspiration from Wang and Bertsekas,  (2016) and further subsample a fixed
number of constraints uniformly at random at each step of the optimization.
For convex problems, this form of constraint subsampling can be shown to
converge to the optimal solution (Wang and Bertsekas, , 2016). Though global
guarantees are elusive for non-convex problems, we follow a long tradition of
machine learning research that achieves impressive results on complex tasks
with stochastic gradient descent. We record the largest (maximization) and
smallest (minimization) values over multiple runs with different random
initializations.

##  4 Experiments

Since ground truth causal effects must be known to properly evaluate the
validity of our bounds, we primarily make use of simulated datasets. All
experiments in the main paper are in a partially identifiable setting. As a
sanity check, we show experiments for identifiable settings as well as for
other datasets and constraint formulations in Appendix F.

###  4.1 Treatment choices

When visualizing results for multidimensional treatments $X$, we vary the
interventional values $x^{\star}$ along a single treatment dimension, keeping
the remaining components at fixed values. While this allows us to show
continuous treatment effect curves, we note our method can compute bounds for
any multidimensional intervention $do(X=x^{\star})$. Specifically, in the
results shown here, we vary the first component of $X$ and fix the values of
the other components to their empirical marginal means. For each figure, we
include a kernel density estimate of the marginal distribution of the
empirically observed treatments to distinguish “data poor” from “data rich”
regions. In “data poor” regions (towards the tails of the empirical
distribution), we may expect our bounds to become looser, as less information
about the data-matching constraints is available.

###  4.2 Baselines

We compare our method to a naïve regression and the method of Hu et al.,
(2021). We are unaware of any other competitor capable of estimating bounds in
the presence of multivariate, continuous treatments/instruments/mediators.

  * •

Regression with MLP. We naïvely fit an MLP with quadratic loss to predict
outcome $Y$ from the multidimensional treatment $X$, a modeling approach that
assumes no confounding at all.

  * •

GAN framework. Hu et al.,  (2021) parameterize the causal model (i.e., its
exogenous random variables and structural equations) with neural networks and
apply the adversarial learning framework to search the parameter space. We
adjust the model in their code to our examples, but otherwise leave
hyperparameter choices and convergence criteria untouched. The error allowed
on matching the observed distribution is $\eta=0.001$ (not our $\eta$) and
bounds are calculated as the mean over the effects of the 50 last models with
error below $\eta$, plus/minus the variance of those 50 values.

###  4.3 Implementation choices

Choice of response functions. As described in Equation 4, we choose linear
combinations of non-linear basis functions $\\{\psi_{k}\\}_{k\in[K]}$. For our
experiments we mostly work with a set of $K$ _neural basis functions_ ,
obtained from the last hidden-layer activations of an MLP fit to the observed
data $\\{(x_{i},y_{i})\\}_{i\in[n]}$, as well as (multivariate) polynomials up
to a fixed degree (see Section E.4 for details). For two-dimensional
treatments, we use $K=6$ and $K=3$ for polynomial and neural basis functions,
respectively. For three-dimensional treatments we use $K=10$ for both
polynomial and neural basis functions. This allows up to quadratic terms in
the polynomial basis for both two- and three-dimensional treatments.

For our method, we individually compute lower ($\bigtriangledown$) and upper
($\triangle$) bounds at multiple values $x^{\star}\in\mathbb{R}$ for one
dimension of the treatment and show how it compares to the true causal effect
( ) and naïve regression ( ). The lines shown for lower and upper ( ) bounds
are univariate cubic splines fit to the bounds for individual $x^{\star}$-grid
values. We use $n=10,000$ i.i.d. sampled datapoints for each experiment and
subsample $M=100$ datapoints uniformly at random for the data constraints
[c-data]. For the $X\,|\,Z$ model in the IV setting we use a quadratic
conditional spline normalizing flow. Details of the implementation, including
hyperparameters, are described in Appendix E.

Getting final bound estimates. For each $x^{\star}$, we run 5 optimizations
with different seeds each for the lower and upper bounds. For the final upper
bound estimate, we take the maximum of the 5 upper bound estimates we get, and
for the final lower bound estimate we take the minimum of the 5 lower bound
estimates we get. We follow the same process for the GAN bounds.

###  4.4 Results

We focus on datasets where the true effect is a linear or quadratic polynomial
function of the treatment $X$. A full glossary of datasets can be found in
Appendix D, along with exact structural equations for each setting. While the
naming of the datasets is intuitive, Figure 4 provides a short description of
the naming logic. We now describe the results obtained for selected datasets
with more results and details in Appendix F.

IV-lin-2d-strong. This dataset is simulated from an IV setting, where the true
effect $Y\,|\,do(X)$ is linear in a two-dimensional treatment $X$. The naïve
regression $\operatorname{\mathbb{E}}[Y\,|\,X]$ differs substantially from the
true effect, indicating strong confounding (see, Figure 2 (bottom left)).

We start by assessing the effect of the flexibility of response functions on
our bounds in Figure 2 (bottom left). Choosing less flexible basis functions
(polynomials) yields tighter bounds, highlighting the flexibility of our
approach in obtaining more informative bounds when more restrictive
assumptions are made. Our method can also accommodate alternative constraint
formulations and slack parameters. We show experiments demonstrating this
flexibility in Appendix F. While sometimes loose, especially for flexible
basis functions (MLP) and in “data poor” regions towards the tails of the
empirical distribution of observed treatments, our bounds contain the true
causal effect for all $x^{\star}$. This validity holds for all other
experiments as well.

Next, Figure 2 (bottom right) compares our approach to the GAN framework
baseline. First, we note that methods that fit a model and bounds at the same
time, like Hu et al.,  (2021), essentially need to solve a bilevel
optimization problem, where the same parameters are shared in the definition
of the constraints (which require optimizing a measure of fitness) and the
causal quantity of interest (which requires optimizing a different problem).
Such a strategy can be wasteful and potentially more prone to optimization and
regularization challenges. Our goal in comparing333For generality, we use the
_neural basis functions_ for all experiments where we compare to the GAN
framework, as they are expected to give wider bounds than polynomials. to Hu
et al.,  (2021) is to show that our framework simplifies constraints through
structural encoding and avoids full density estimation while giving comparable
or better results than those from the GAN framework. Indeed, Figure 2 (bottom
right) shows that despite the flexible neural basis functions, our approach
can yield tighter bounds and avoid the instabilities observed in the GAN
approach when we move towards the tails of the observational distribution.

We show similar comparisons for linear and quadratic true effect IV settings
with strong and weak confounding as well as two- and three-dimensional
treatments in Appendix F.

IV-lin-2d-weak. This is an IV setting with weak confounding. In Figure 2 (top
left) we show how our bounds behave under different choices of
$\operatorname{\mathrm{dist}}$, namely under the entry-wise
$\|\cdot\|_{\infty,\infty}$ norm (which results in $M\cdot L$ constraints in
the augmented Lagrangian) versus the entry-wise $\|\cdot\|_{2,2}$ norm
(yielding only a single constraint). The obtained bounds are compatible and
comparable, indicating relatively mild effects of the choice of the data-
matching criterion $\operatorname{\mathrm{dist}}$.

LM-lin-2d-strong. Figure 2 (top right) shows results for data from a leaky
mediator setting, where the true effect is a linear function of the mediator
$M$, and the treatment $X$ is two-dimensional. Confounding between the
mediator and the effect is relatively strong. These results corroborate our
findings regarding stability and reliability of our bounds compared to the GAN
framework from the IV setting, despite using the neural basis function. Again,
additional results for different structural assumptions, treatment dimensions,
and confounding strength are presented in Appendix F.

##  5 Discussion

Causal modeling inevitably involves a trade-off between the strength of input
assumptions and the specificity of resulting inferences. If assumptions do not
correspond to reality, the inference that follows is unwarranted. The role of
a method is to provide the sandbox in which assumptions can be transparently
expressed and informative constraints, where available, are not wasted.
Ultimately, no automated method can choose the untestable assumptions on which
any particular application of causal inference relies. See Appendix A for
further discussion.

After parameterization, the optimization problem in Equation 1 will generally
result in a _non-convex objective with non-convex constraints_. Therefore, our
proposed gradient-based local optimization may not converge to a global
optimum, possibly rendering our bounds overly tight. Empirically, we do not
observe evidence of consistently getting stuck in bad local optima. Because we
optimize both bounds individually for each value of $x^{\star}$, our method
may be computationally expensive when the intervention space is very large.
However, it is well-suited for scenarios where we want reliable bounds on a
well-defined set of plausible interventions. Finally, we have not accounted
for the uncertainty of our bounds. Confidence or credible intervals for both
extrema can help practitioners evaluate the reliability of causal inferences,
and are an interesting direction for future work.

From an ethical perspective, more reliable methods for bounding causal effects
could help quantify causal effects of high-stakes or consequential decisions
when identifiability is impossible. They may also alleviate the temptation to
blindly accept whichever assumptions are necessary to ensure point
identification. However, by the same token, bounds may still misrepresent
treatment effects when untestable assumptions are not satisfied. We provide
more details on ethical considerations in Appendix H.

We have introduced a stochastic causal program for bounding treatment effects
in partially identifiable settings. Our approach does not rely on the typical
assumptions of linearity, monotonocity, or additivity. We presented a novel
parameterization for decoupling observed from unobserved potential outcomes,
and derived an efficient procedure for optimizing over a set of causal models
consistent with the data. Experiments demonstrate that our method produces
valid and informative bounds in a wide range of settings, including with
continuous and multidimensional instruments, mediators, treatments, and
outcomes.

## References

  * Angrist and Imbens, (1995) Angrist, J. D. and Imbens, G. W. (1995).  Two-stage least squares estimation of average causal effects in models with variable treatment intensity.  J. Am. Stat. Assoc., 90(430):431–442.
  * Balke and Pearl, (1994) Balke, A. and Pearl, J. (1994).  Counterfactual probabilities: Computational methods, bounds and applications.  In Proceedings of the 10th Conference on Uncertainty in Artificial Intelligence, pages 46–54.
  * Balke and Pearl, (1997) Balke, A. and Pearl, J. (1997).  Bounds on treatment effects from studies with imperfect compliance.  J. Am. Stat. Assoc., 92(439):1171–1176.
  * Bennett et al., (2019) Bennett, A., Kallus, N., and Schnabel, T. (2019).  Deep generalized method of moments for instrumental variable analysis.  In Advances in Neural Information Processing Systems.
  * Bica et al., (2021) Bica, I., Alaa, A. M., Lambert, C., and van der Schaar, M. (2021).  From real-world patient data to individualized treatment effects using machine learning: Current and future methods to address underlying challenges.  Clin. Pharmacol. Ther., 109(1).
  * Bingham et al., (2019) Bingham, E., Chen, J. P., Jankowiak, M., Obermeyer, F., Pradhan, N., Karaletsos, T., Singh, R., Szerlip, P. A., Horsfall, P., and Goodman, N. D. (2019).  Pyro: Deep universal probabilistic programming.  J. Mach. Learn. Res., 20:28:1–28:6.
  * Chernozhukov et al., (2018) Chernozhukov, V., Chetverikov, D., and Kato, K. (2018).  Inference on causal and structural parameters using many moment inequalities.  Rev. Econ. Stud., 86(5):1867–1900.
  * Chickering and Pearl, (1996) Chickering, D. M. and Pearl, J. (1996).  A clinician’s tool for analyzing non-compliance.  In Proceedings of the 13th AAAI Conference on Artificial Intelligence, pages 1269–1276.
  * Daubechies et al., (2021) Daubechies, I., DeVore, R., Foucart, S., Hanin, B., and Petrova, G. (2021).  Nonlinear approximation and (deep) relu networks.  Constr. Approx., pages 1–46.
  * Dawid, (2003) Dawid, A. (2003).  Causal inference using influence diagrams: the problem of partial compliance.  In Green, P., Hjort, N., and Richardson, S., editors, Highly Structured Stochastic Systems, pages 45–65. Oxford University Press.
  * Drton and Richardson, (2008) Drton, M. and Richardson, T. (2008).  Graphical methods for efficient likelihood inference in Gaussian covariance models.  J. Mach. Learn. Res., 9:893–914.
  * Duarte et al., (2021) Duarte, G., Finkelstein, N., Knox, D., Mummolo, J., and Shpitser, I. (2021).  An automated approach to causal inference in discrete settings.  _arXiv_ preprint, 2109.13471.
  * Gentzel et al., (2021) Gentzel, A. M., Pruthi, P., and Jensen, D. (2021).  How and why to use experimental data to evaluate methods for observational causal inference.  Proceedings of the 38th International Conference on Machine Learning, pages 3660–3671.
  * Gunsilius, (2019) Gunsilius, F. (2019).  Bounds in continuous instrumental variable models.  _arXiv_ preprint, 1910.09502.
  * Gunsilius, (2021) Gunsilius, F. F. (2021).  Nontestability of instrument validity under continuous treatments.  Biometrika, 108(4):989–995.
  * Gupta et al., (2020) Gupta, M., Louidor, E., Mangylov, O., Morioka, N., Narayan, T., and Zhao, S. (2020).  Multidimensional shape constraints.  Proceedings of the 37th International Conference on Machine Learning, pages 3918–3928.
  * Hartford et al., (2017) Hartford, J., Lewis, G., Leyton-Brown, K., and Taddy, M. (2017).  Deep IV: A flexible approach for counterfactual prediction.  In Proceedings of the 34th International Conference on Machine Learning.
  * Hernán and Robins, (2020) Hernán, M. A. and Robins, J. M. (2020).  Causal Inference: What If.  Chapman & Hall/CRC Studies, Boca Raton.
  * Hu et al., (2021) Hu, Y., Wu, Y., Zhang, L., and Wu, X. (2021).  A generative adversarial framework for bounding confounded causal effects.  In Proceedings of the 35th AAAI Conference on Artificial Intelligence.
  * Imbens and Rubin, (2015) Imbens, G. W. and Rubin, D. B. (2015).  Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction.  Cambridge University Press, Cambridge.
  * Kallus and Zhou, (2020) Kallus, N. and Zhou, A. (2020).  Confounding-robust policy evaluation in infinite-horizon reinforcement learning.  In Advances in Neural Information Processing Systems.
  * Kilbertus et al., (2019) Kilbertus, N., Ball, P. J., Kusner, M. J., Weller, A., and Silva, R. (2019).  The sensitivity of counterfactual fairness to unmeasured confounding.  In Proceedings of the 35th Conference on Uncertainty in Artificial Intelligence.
  * Kilbertus et al., (2020) Kilbertus, N., Kusner, M. J., and Silva, R. (2020).  A class of algorithms for general instrumental variable models.  In Advances in Neural Information Processing Systems.
  * Kilbertus et al., (2017) Kilbertus, N., Rojas Carulla, M., Parascandolo, G., Hardt, M., Janzing, D., and Schölkopf, B. (2017).  Avoiding discrimination through causal reasoning.  Advances in neural information processing systems, 30.
  * Kusner et al., (2017) Kusner, M. J., Loftus, J., Russell, C., and Silva, R. (2017).  Counterfactual fairness.  Advances in neural information processing systems, 30.
  * Manski, (1990) Manski, C. F. (1990).  Nonparametric bounds on treatment effects.  Am. Econ. Rev., 80(2):319–323.
  * Muandet et al., (2020) Muandet, K., Mehrjou, A., Lee, S. K., and Raj, A. (2020).  Dual instrumental variable regression.  In Advances in Neural Information Processing Systems, volume 33, pages 2710–2721.
  * Nocedal and Wright, (2006) Nocedal, J. and Wright, S. (2006).  Numerical optimization.  Springer Science & Business Media.
  * Papamakarios et al., (2021) Papamakarios, G., Nalisnick, E., Rezende, D. J., Mohamed, S., and Lakshminarayanan, B. (2021).  Normalizing flows for probabilistic modeling and inference.  J. Mach. Learn. Res., 22(57):1–64.
  * Paszke et al., (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S. (2019).  Pytorch: An imperative style, high-performance deep learning library.  In Advances in Neural Information Processing Systems 32, pages 8024–8035.
  * Pearl, (2009) Pearl, J. (2009).  Causality.  Cambridge University Press, New York.
  * Richardson, (2003) Richardson, T. (2003).  Markov properties for acyclic directed mixed graphs.  Scand. J. Stat., 30(1):145–157.
  * Robins et al., (2015) Robins, J. M., VanderWeele, T. J., and Gill, R. D. (2015).  A proof of Bell’s inequality in quantum mechanics using causal interactions.  Scandinavia Journal of Statistics, 42:329–335.
  * Schölkopf et al., (2022) Schölkopf, B., Herbrich, R., and Smola, A. J. (2022).  A generalized representer theorem.  Lecture Notes in Computer Science, 2111:416–426.
  * Shen et al., (2021) Shen, Z., Yang, H., and Zhang, S. (2021).  Neural network approximation: Three hidden layers are enough.  Neural Netw., 141:160–173.
  * Singh et al., (2019) Singh, R., Sahani, M., and Gretton, A. (2019).  Kernel instrumental variable regression.  In Advances in Neural Information Processing Systems, pages 4595–4607.
  * Sonja A. Swanson et al., (2018) Sonja A. Swanson, M. A. H., Miller, M., Robins, J. M., and Richardson, T. S. (2018).  Partial identification of the average treatment effect using instrumental variables: Review of methods for binary instruments, treatments, and outcomes.  J. Am. Stat. Assoc., 113:933–947.
  * Spirtes et al., (2000) Spirtes, P., Glymour, C., and Scheines, R. (2000).  Causation, Prediction and Search.  MIT Press.
  * Tchetgen et al., (2022) Tchetgen, E. J. T., Ying, A., Cui, Y., Shi, X., and Miao, W. (2022).  An introduction to proximal causal learning.  arXiv:2009.10982.
  * Vandenbroucke et al., (2016) Vandenbroucke, J. P., Broadbent, A., and Pearce, N. (2016).  Causality and causal inference in epidemiology: the need for a pluralistic approach.  Int. J. Epidemiol., 45(6):1776–1786.
  * Wang and Bertsekas, (2016) Wang, M. and Bertsekas, D. P. (2016).  Stochastic first-order methods with random constraint projection.  SIAM J. Optim., 26(1):681–717.
  * Wasserman, (2006) Wasserman, L. (2006).  All of Nonparametric Statistics.  Springer.
  * Wolfe et al., (2019) Wolfe, E., Spekkens, R. W., and Fritz, T. (2019).  The inflation technique for causal inference with latent variables.  J. Causal Inference, 7(2).
  * (44) Wu, Y., Zhang, L., and Wu, X. (2019a).  Counterfactual fairness: Unidentification, bound and algorithm.  In Proceedings of the 28th International Joint Conference on Artificial Intelligence, pages 1438–1444.
  * (45) Wu, Y., Zhang, L., Wu, X., and Tong, H. (2019b).  Pc-fairness: A unified framework for measuring causality-based fairness.  In Advances in Neural Information Processing Systems.
  * Xia et al., (2021) Xia, K., Lee, K.-Z., Bengio, Y., and Bareinboim, E. (2021).  The neural-causal connection: expressiveness, learnability, and inference.  In Advances in Neural Information Processing Systems.
  * Zhang and Bareinboim, (2020) Zhang, J. and Bareinboim, E. (2020).  Designing optimal dynamic treatment regimes: A causal reinforcement learning approach.  In Proceedings of the 37th International Conference on Machine Learning.
  * Zhang and Bareinboim, (2021) Zhang, J. and Bareinboim, E. (2021).  Bounding causal effects on continuous outcomes.  In Proceedings of the 35th AAAI Conference on Artificial Intelligence.
  * Zhang et al., (2021) Zhang, J., Tian, J., and Bareinboim, E. (2021).  Partial counterfactual identification from observational and experimental data.  _arXiv_ preprint, 2110.05690.
  * Zheng et al., (2018) Zheng, X., Aragam, B., Ravikumar, P. K., and Xing, E. P. (2018).  Dags with no tears: Continuous optimization for structure learning.  In Advances in Neural Information Processing Systems 31.

##  Appendix A On the Choice of Function Space

The function space in Eq. (4) is, at a high-level view, a standard finite
basis expansion with or without adaptive bases. We chose this particular
function space template as it is (i) general, (ii) computationally
advantageous and (iii) scientifically advantageous.

It is general, as it can arbitrarily approximate any practical function of
interest (keeping in mind that allowing for an unbounded number of
discontinuities will result in vacuous bounds (Gunsilius, , 2021)),
particularly when we are dealing with a small number of causal parents. As a
matter of fact, (Gunsilius, , 2019) uses as basis functions a truncated
wavelet expansion. We could likewise use the dictionary of regression spline
basis functions, with splits decided by any rule of interest, including
uniform grids or those based on the training points.

As a matter of fact, no computable estimator of nonparametric functions is
infinite dimensional at any given sample size. In case of kernel machines, for
instance, the celebrated Representer Theorem of (Schölkopf et al., , 2022)
shows how any resulting estimate will be a function that is parameterized
explicitly by the data (and such will be the case for the broad class of
linear smoothers (Wasserman, , 2006, Ch. 5)). As our causal models are _not_
fit to the data directly, but to estimated functionals of the observed
distribution, which we use as constraints in a mathematical program, we can
directly parameterize it in terms of a finite basis expansion.

Moreover, there is nothing stopping us from making each element $\psi_{i}$ of
the basis functions parameterized by, say, a linear transformation of the
function inputs followed by some non-linearity (effectively, a hidden unit in
a feedforward neural network). An entire neural network could also be used for
each $\psi_{i}$, so that the function space becomes a deep model, with the
distribution over the function space confined only to randomness on the
$\theta$ parameters. Nothing in the main algorithm changes, except for extra
terms added to the overall gradient.

This leads to the fact that our function space is computationally
advantageous, particularly when we can eliminate bidirected edges using ideas
from Drton and Richardson,  (2008), combined with the decision to use only the
first two moments of target causal factors. In this case, all constraints can
be computed in closed-form and we can perform a relatively simple Monte Carlo
approximation for the objective function. As a matter of fact, we did not have
to make any distributional assumptions about $\theta$ at all except for the
existence of its first and second moments. There is much to appreciate in the
generic blackbox of Hu et al.,  (2021), but solving full density estimation
problems and requiring a complex nesting of multiple optimization and
simulation steps may not be the preferred way to go in some real-world
applications.

Moreover, response functions only need to be structural with respect to its
manipulable inputs. If we want to model, say, $\mathbb{E}[Y~{}|~{}do(x),w]$
for some pre-treatment covariates $W$, there is no need at all to model $W$
structurally. The causal model can be defined in terms of the response
function $f_{(w,u_{y})}(x)$, where $W$ is just another source of background
variability (though observable, unlike $U_{Y}$). This means that the
_distribution_ of $\theta$ will be parameterized as a black-box function that
includes $w$ as inputs, while the basis function expansion can remain
flexible, and informative, in the small dimensional space of $X$.

Finally, our setup is scientifically advantageous. A naïve view of causal
inference postulates that since functional constraints are untestable
assumptions that will imply invalid bounds if misspecified, then we must, on
purpose, impoverish the language in which we allow causal models to be
described by removing this choice. This line of reasoning is incoherent on
many accounts. First, practitioners do want to express functional constraints,
as partial identification is notoriously uninformative in many cases. Shape
constraints may be fruitful even in purely predictive tasks (e.g Gupta et al.,
, 2020), as further background knowledge to reduce degrees of freedom that
cannot be fully decided by data alone. Many domains have no clear reason to
allow for dose-response curves that have many inflection points, and in fact
the community has been criticized by creating benchmarks with artificial
outcome functions for the sake of illustrating overly complex average
treatment effect estimators, which hardly deliver any advantages in several
empirical domains (e.g., Gentzel et al., , 2021). The argument that
restricting function spaces may lead to invalid bounds is vacuous and misses
the point of causal modeling. A model is what we remove from the space of all
possibilities, that is, the constraints we adopt.444To use a quote usually
apocryphally attributed to Michelangelo: “It is the sculptor’s power, so often
alluded to, of finding the perfect form and features of a goddess, in the
shapeless block of marble; and his ability to chip off all extraneous matter,
and let the divine excellence stand forth for itself.” _The Methodist
Quarterly Review_ (1858), emphasis added. The missing edges in a causal graph
are just that: constraints. Invalid bounds are a consequence of an
overconstrained feasible region, and the function space is part of the game as
much as any decisions about missing edges. The promotion of a type of
constraint (independence constraints) as having a protected place in causal
modeling should not be taken seriously by methodologists or practitioners.

Second, like missing edges, function spaces can be partially tested.
Testability is predictated on assumptions that themselves may be untestable:
in the case of purely missing edges, faithfulness (Spirtes et al., , 2000) and
its variants; in the case of function spaces and the independence constraints
that go along with them, we can test whether our mathematical program has a
feasible solution or not. Passing the mathematical program is not a guarantee
of correctness – which is not surprising, given that e.g. hypothesis testing
and many frameworks of scientific falsifiability are about deciding not to
reject a postulated model, not about confirming it. As another example,
passing Fisher’s sharp null hypothesis by itself does not prove that
individual treatment effects are exactly, or even close to, zero (or whether
counterfactuals have a physical meaning at all!), and so on. Nevertheless,
there are important applications of falsifying a causal model with partial
identifiability (e.g Wolfe et al., , 2019; Robins et al., , 2015), which are
useful so long as it is understood that “testability” is always an asymmetric
concept.

Among the practical uses of constrained function spaces, for instance, is the
possibility of expressing knowledge as smooth deviances from responses that do
not control for unmeasured confounding: e.g., parameterizing
$Y=f_{(w,u_{y})}(x)$ as $Y=\mathbb{E}[Y~{}|~{}w,x]+r_{(w,u_{y})}(x)$, the
latter function $r_{(w,u_{y})}(x)$ taking particularly smooth shapes.

One line of criticism that deserves more serious consideration was raised by
Dawid Dawid,  (2003) in the context of classical categorical models for
instrumental variables: postulating constraints on structural
equations/response functions is equivalent to making counterfactual claims
which are never directly testable, even with perfect randomized controlled
trials (RCTs). This motivated Dawid to create an alternative to partial
identification in discrete IV settings without postulating the existence of
response functions. A discussion of interventional (response-function-free)
and counterfactual (adopting response functions) models is provided by Sonja
A. Swanson et al.,  (2018). In our case, expressing models in terms of latent
variables and conditional distributions falls back to the issue we discussed
in the main text, by which it is unclear how to interpret an infinite
dimensional latent variable space and its relation to the conditional
distributions, and all the computational complications that follow. Our take
is that the response functions are idealizations that do not need to
correspond to counterfactuals if all that is required is how well we can model
interventional distributions – themselves testable implications in the sense
of comparisons against RCTS. Without RCTs, the formulation has a degree of
falsifiability in terms of the feasible region of the mathematical program, as
discussed above.

##  Appendix B Leaky mediator

In the “leaky mediator” setup, we have that the effect of treatment $X$ on
outcome $Y$ is completely mediated by some vector $M$. On top of that, $X$ and
$Y$ are confounded, and so are $M$ and $Y$. Without the latter, this causal
structure would satisfy the front-door criterion, and be nonparametrically
identifiable via the do-calculus (Pearl, , 2009). A causal diagram with our
proposed parameterization is shown in Figure 1(d), right column.

For fixed background variables $U_{Y}$ encoded indirectly as function
parameters $\theta_{y}$, the structural equation for $Y$ is given by
$f_{\theta_{y}}(m)$. Assume we represent the density function of $M$ given
$X$, $p_{m}(\cdot\,|\,x)$, with an invertible conditional normalizing flow
$M=h_{X}(N_{M})$, where $N_{M}\sim\mathcal{N}(0,1)$. In what follows, we adopt
no model for $p_{x}(x)(=p(n_{x}))$, using its empirical distribution
$\hat{p}_{x}(\cdot)$ as an estimate.

Representing the objective function. A main difference between this and the IV
model is that we also need to marginalize the mediator $M$. We can interpret
the expectation $\operatorname{\mathbb{E}}[Y\,|\,do(x^{\star})]$ as following
the same generative model of the observational distribution, but where node
$X$ gets replaced by constant $x^{\star}$, with edge $N_{X}\rightarrow X$ also
removed. This leads to the following integral and its Monte Carlo
approximation:

$\displaystyle\mathbb{E}[Y\,|\,do(x^{\star})]$ $\displaystyle=\int
f_{\theta_{y}}(m)p(m,n_{x},\theta_{y}\,|\,do(x^{\star}))\,dm\,dn_{x}\,d\theta_{y}$
$\displaystyle\approx\frac{1}{n}\sum_{i}\theta_{y}^{(i)\top}\psi(m^{(i)}),$
(10)

where

$\begin{array}[]{rcl}n_{x}^{(i)}\,|\,do(x^{\star})&\sim&\hat{p}_{x}(\cdot)\\\
m^{(i)}\,|\,n_{x}^{(i)},do(x^{\star})&\sim&\hat{p}_{m}(\cdot\,|\,x^{\star})\\\
n_{m}^{(i)}\,|\,m^{(i)},do(x^{\star})&=&h_{x^{\star}}^{-1}(m^{(i)})\\\
\theta_{y}^{(i)}\,|\,n_{m}^{(i)},m^{(i)},n_{x}^{(i)},do(x^{\star})&\sim&p_{\eta}(\cdot\,|\,n_{x}^{(i)},n_{m}^{(i)})\end{array}$

Here, $n$ is the number of Monte Carlo samples; $\hat{p}_{x}(\cdot)$ is the
empirical distribution of $X$ (and, hence, sampling from it is just sampling
with replacement from the training set for $X$). This can be interpreted
conceptually, with the actual algorithm being the sampling of $N_{M}$ from a
standard Gaussian followed by the computation of $M$ from $x^{\star}$ and
$N_{M}$. As in the main text, $p_{\eta}(\cdot~{}|~{}\cdot,\cdot)$ only needs
to be specified up to its first and second moments.

Estimating the constraints. Our optimization is over multivariate Gaussian
distributions of $\theta_{y}$ with the mean and covariance modeled by MLPs
with parameters $\eta_{0}$ and $\eta_{1}$, respectively. Let’s say the MLPs
are $\mu_{\eta_{0}}(n_{x},n_{m})$ and $\Sigma_{\eta_{1}}(n_{x},n_{m})$. The
constraints of the optimization are to match the generated data to the
observed distribution, for example by matching the first and second moments of
$Y\,|\,X,M$ as in the IV case.

$\displaystyle\mathbb{E}[Y\,|\,x,m]$ $\displaystyle=\int
f_{\theta_{y}}(m)p(\theta_{y}\,|\,x,m)\,d\theta_{y}$ (11)
$\displaystyle=\mathbb{E}[\theta^{\top}\,|\,x,m]\psi(m)$ (12)
$\displaystyle=\mu_{\eta_{0}}(x,n_{m}(x,m))^{\top}\psi(m)$ (13)

This is analogous to the IV setting because of the similar factorization of
the background variables as a DAG after re-encoding the bidirected edges as
directed edges. After transforming the causal background variables $U_{X}$,
$U_{M}$ into probabilistic “seeds” $N_{X}$, $N_{M}$ for $p(x)$ and
$p(m\,|\,x)$ , this results in the factorization
$p(n_{x})\,p(n_{m})\,p(\theta_{y}\,|\,n_{m},n_{x})$. As $N_{X}$ is identically
distributed as $X$ and so can be estimated separately, and $p(n_{m})$ is fixed
by design, the same general pipeline follows.

##  Appendix C Extensions

We introduced our problem formulation as a very generic definition of model
proximity in terms of some $\operatorname{\mathrm{dist}}(p_{\eta},\hat{p})$,
and in particular the formulation

$\operatorname{\mathrm{dist}}(p_{\eta},\hat{p})\equiv||f^{p_{\eta}}-f^{\hat{p}}||_{\infty},$

where $f^{\hat{p}}$ is a finite-dimensional vector of functionals of the
(estimated) observable distribution, and $f^{p_{\eta}}$ is a finite-
dimensional vector of functionals of the distribution as implied by an
unidentifiable structural causal model. Equation 6 is a particular
implementation of this idea. Due to the challenge of evaluating this norm, we
suggest Monte Carlo approximation for the optimization procedure, based on
Wang and Bertsekas,  (2016).

We have also presented results for alternative distance metrics, including
special cases such as metrics that can be associated with a single Lagrange
multiplier such as $||f^{p_{\eta}}-f^{\hat{p}}||_{2}$. We note that also
combinations of a limited number of metrics are possible such as

$\max\\{||f^{p_{\eta}}_{1}-f^{\hat{p}}_{1}||_{1},||f^{p_{\eta}}_{1}-f^{\hat{p}}_{2}||_{1},\dots,||f^{p_{\eta}}_{K}-f^{\hat{p}}_{K}||_{1}\\},$

where each $k=1,2,\dots,K$ describes a cluster of functionals that are easy to
compare on the same scale.

We emphasize an approach where $f^{\hat{p}}$ is learned prior to solving an
optimization problem. It may be of interest to use the raw data directly and
try to fit $p_{\eta}$ while searching, among all equally best fitting models,
for the one with the highest/lowest corresponding causal effect. There are
some appealing properties of using something like (Monte Carlo approximations
of) the log-likelihood function evaluated at the training points (i.e., the
negative KL divergence between model and data), maximum mean discrepancy
measures, or generative adversarial network criteria, as they are global
measures of goodness-of-fit. Future work on such approaches should also be
considered. However, it is not entirely obvious how nested optimization would
complicate the solution to the problem, as optimizing (and regularizing!) the
fit plays the role of a constraint for the “true” objective function that is
optimized (the causal effect). Moreover, fitting a generative model for the
entire joint distribution is an obviously undesirable nuisance task if we can
get essentially the same bounds by estimating only one or two moments of
interest.

Another extension is the provision of a calculus for a minimal
parameterization of arbitrary structural causal models for a given causal
query of interest. Looking back at our instrumental variable scenario, we
could have created a single set of latent variables $N_{XY}$ that would be
common causes of $X$ and $Y$ and model the structural equations for both $X$
and $Y$ at the same time, instead of fixing $p(X\,|\,Z)$ a priori and
deterministically extracting latent variables $N$ from $X$ and $Z$. In
general, we could add independent latent variables for each clique in the
bidirected graph component of the causal graph. This, however, is very
wasteful. There is no need to create a causal model for the $(Z,X)$ marginal –
which not only would assume that $Z$ is an unconfounded cause of $X$, which is
unnecessary, but would also waste computation and stability trying to match
the observable marginal of $(Z,X)$ to the corresponding causal-model-implied
marginal, which is also completely unnecessary. Given that $Z$ and $X$ can be
high-dimensional while $Y$ is a scalar, this is clearly a bad idea. Therefore,
the structure-blind strategy of creating latent variables for each bidirected
clique is convenient but not ideal, and a smarter automated way of generating
minimal causal parameterizations for arbitrary graphs and causal queries is
needed.

##  Appendix D Glossary of Synthetic Datasets

We describe here the synthetic datasets we use in our experiments.

Instrumental variable Leaky mediation $Z$$\in\mathbb{R}^{p}$
instruments$X$$\in\mathbb{R}^{q}$ treatments$Y$$\in\mathbb{R}$ outcome$C$
$X$$\in\mathbb{R}^{p}$ treatments$M$$\in\mathbb{R}^{q}$
mediators$Y$$\in\mathbb{R}$ outcome$U$$C$

Figure 3: These are the structural equations we work with. In the description
of the datasets, the dimensions of each variable are denoted by a subscript.
For instance if $X$ is $2$-dimensional, we write it as $X=(x_{1},x_{2})$.

IV

\- lin \- 2d \- strong \- add (optional)

Figure 4: Naming logic for the datasets. The first segment mentions whether it
is the Instrumental Variable (IV) or Leaky Mediator (LM) setting. The second
segment says whether $y$ is linear or quadratic in $x$. The third segment
tells the dimension of the treatment. The fourth segment denotes the strength
of the confounding, strong or weak. The fifth and last segment is optional,
and says whether the confounding is additive, which in the IV setting makes
the effect identifiable.

We use various polynomial datasets where the causal effect is a polynomial
function of a single or multidimensional treatment. We provide the
construction of each of these here. Figure 3 shows the structural graphs. In
both graphs, every node except $Y$ (outcome) could be multi-dimensional. If a
node, say $C$, is multi-dimensional, we index the dimension with a subscript.
That is, we write $C=(c_{1},c_{2})$. Figure 4 visually describes the naming
logic behind the datasets to make the exposition clearer.

###  D.1 Scalar treatment

We now describe settings where the treatment $X$ is a scalar,which we present
as a sanity check. The noises, confounder, and instrument follow

$e_{x},e_{y},c,z\sim\mathcal{N}(0,1)$

  * •

IV-lin-1d-weak-add ($f_{y}$ linear in $x$, weak additive confounding)

$\displaystyle f_{x}(z,c,e_{x})$ $\displaystyle=3z+0.5c+e_{x}$ $\displaystyle
f_{y}(x,c,e_{y})$ $\displaystyle=x-6c+e_{y}$

  * •

IV-quad-1d-strong ($f_{y}$ quadratic in $x$, strong non-additive confounding)

$\displaystyle f_{x}(z,c,e_{x})$ $\displaystyle=0.5z+3c+e_{x}$ $\displaystyle
f_{y}(x,c,e_{y})$ $\displaystyle=0.3x^{2}-1.5xc+e_{y}$

  * •

IV-quad-1d-weak ($f_{y}$ quadratic in $x$, weak non-additive confounding)

$\displaystyle f_{x}(z,c,e_{x})$ $\displaystyle=3z+0.5c+e_{x}$ $\displaystyle
f_{y}(x,c,e_{y})$ $\displaystyle=0.3x^{2}-1.5xc+e_{y}$

###  D.2 IV model

We now describe datasets satisfying the IV assumptions.

2D treatment. The noises, confounder, and instruments follow

$\displaystyle c,z,e_{x}$ $\displaystyle\sim\mathcal{N}^{2}(0,1)$
$\displaystyle e_{y}$ $\displaystyle\sim\mathcal{N}(0,1)$

  * •

IV-lin-2d-strong ($f_{y}$ linear in $x$, strong non-additive confounding)

$\displaystyle f_{x}(z,c,e_{x})$ $\displaystyle=0.5z+2c+e_{x}$ $\displaystyle
f_{y}(x,c,e_{y})$
$\displaystyle=x_{1}+x_{2}-3(x_{1}+x_{2})(c_{1}+c_{2})+e_{y}$

  * •

IV-lin-2d-weak ($f_{y}$ linear in $x$, weak non-additive confounding)

$\displaystyle f_{x}(z,c,e_{x})$ $\displaystyle=2z+c+e_{x}$ $\displaystyle
f_{y}(x,c,e_{y})$ $\displaystyle=5x_{1}+6x_{2}-x_{1}(c_{1}+c_{2})+e_{y}$

  * •

IV-quad-2d-strong-add ($f_{y}$ quadratic in $x$, strong additive confounding)

$\displaystyle f_{x}(z,c,e_{x})$ $\displaystyle=z+2c+e_{x}$ $\displaystyle
f_{y}(x,c,e_{y})$ $\displaystyle=2x_{1}^{2}+2x_{2}^{2}-(c_{1}+c_{2})+e_{y}$

  * •

IV-quad-2d-weak ($f_{y}$ quadratic in $x$, weak non-additive confounding)

$\displaystyle{2}f_{x}(z,c,e_{x})$ $\displaystyle=2z+c+e_{x}$ $\displaystyle
f_{y}(x,c,e_{y})$
$\displaystyle=5x_{1}^{2}+6x_{2}^{2}-(x_{1}+x_{2})(c_{1}+c_{2})+e_{y}$

3D treatment. The noises, confounders, and instruments follow

$\displaystyle c,z,e_{x}$ $\displaystyle\sim\mathcal{N}^{3}(0,1)$
$\displaystyle e_{y}$ $\displaystyle\sim\mathcal{N}(0,1)$

  * •

IV-quad-3d-weak ($f_{y}$ quadratic in $x$, weak non-additive confounding)

$\displaystyle f_{x}(z,c,e_{x})$ $\displaystyle=2z+c+e_{x}$ $\displaystyle
f_{y}(x,c,e_{y})$ $\displaystyle=2x_{1}^{2}+2x_{2}^{2}+2x_{3}$
$\displaystyle\phantom{=}-0.3(x_{2}+x_{3})(c_{1}+c_{2}+c_{3})+e_{y}$

###  D.3 Leaky Mediator

We now describe datasets following the leaky mediator model with noises, and
confounder following

$\displaystyle e_{x},e_{m},c,u$ $\displaystyle\sim\mathcal{N}^{2}(0,1)$
$\displaystyle e_{y}$ $\displaystyle\sim\mathcal{N}(0,1)$

  * •

LM-lin1-2d ($f_{y}$ linear in $x$, strong confounding)

$\displaystyle f_{x}(u,e_{x})$ $\displaystyle=u+e_{x}$ $\displaystyle
f_{m}(x,c,e_{m})$ $\displaystyle=x+3c-e_{m}$ $\displaystyle
f_{y}(m,c,u,e_{y})$ $\displaystyle=2m_{1}+m_{2}$
$\displaystyle\phantom{=}-(m_{1}+m_{2})(c_{1}+c_{2}+u_{1}+u_{2})+e_{y}$

  * •

LM-lin2-2d ($f_{y}$ linear in $x$, weak confounding)

$\displaystyle f_{x}(u,e_{x})$ $\displaystyle=u+e_{x}$ $\displaystyle
f_{m}(x,c,e_{m})$ $\displaystyle=3x+c-e_{m}$ $\displaystyle
f_{y}(m,c,u,e_{y})$ $\displaystyle=2m_{1}+m_{2}-$
$\displaystyle\phantom{=}0.3(m_{1}+m_{2})(c_{1}+c_{2}+u_{1}+u_{2})+e_{y}$

##  Appendix E Implementation Details

We use $n=10\,000$ data points for all our simulations. All implementation is
in Python, using PyTorch (Paszke et al., , 2019). Code for the experiments is
provided with the supplementary material.

###  E.1 Satisfying [c-data]

In Section 3.4 we explained how we match the observed data distribution by
matching scalar statistics for $L$ dictionary functions
$\\{\phi_{l}:\mathcal{Y}\to\mathbb{R}\\}_{l=1}^{L}$. In practice, we take
$L=2$ and match the first and second moments of $Y\,|\,\\{X,Z\\}$, which is to
say that we set $\phi_{1}(Y)=Y$ and $\phi_{2}(Y)=Y^{2}$. We learn $\phi_{1}$
and $\phi_{2}$ by regressing MLPs on the observed data
$\\{x_{i},z_{i},y_{i}\\}_{i\in[n]}$ and
$\\{x_{i},z_{i},y_{i}^{2}\\}_{i\in[n]}$ respectively, with
$\\{x_{i},z_{i}\\}_{i\in[n]}$ as the input and $\\{y_{i}\\}_{i\in[n]}$ and
$\\{y_{i}^{2}\\}_{i\in[n]}$ as the target values. The outputs of these MLPs
are then approximations of $\operatorname{\mathbb{E}}[Y\,|\,X,Z]$ and
$\operatorname{\mathbb{E}}[Y^{2}\,|\,X,Z]$, respectively. A significant
advantage of this approach is that we can evaluate the constraints in closed
form under our construction (see Equation 8). We use MLPs with $3$ hidden
layers of sizes $(64,32,16)$ (first to last) and train it with a batch size of
$512$ for $200$ epochs, with a learning rate of $0.01$. We use ReLU
activations. Identical settings were used for the leaky mediator, with the
only difference being that we match $Y\,|\,\\{M,X\\}$ instead of
$Y\,|\,\\{X,Z\\}$.

GAN comparison

(a) IV-lin-2d-weak.

(b) LM-lin-2d-weak.

Figure 5: Comparison with GAN

###  E.2 Satisfying [c-struct]

We describe in Section 3.3 how we write the mean and variance of $X$ as a
function of $Z$. We considered two distinct ways of modeling $X\,|\,Z$.

Gaussian. . We write the mean and variance of $X$ as a function of $Z$. That
is to say, $X=h_{Z}(N):=A(Z)N+b(Z)$. We parameterize $A(Z)$ as
$A=L^{\top}L+\Omega$. A Cholesky factor $L:\mathcal{Z}\to\mathbb{R}^{p\times
p}$ ensures symmetry and a diagonal matrix $\Omega$ adds small constants to
the diagonal to ensure symmetry and positive definiteness. We then learn
$L^{\top}L+\Omega$ from observed data $\mathcal{D}$ by maximizing the log-
likelihood of $\hat{p}(X\,|\,Z)$. In practice, the parameters of $A(Z)$ as
described above, and the parameters of $b(Z)$ (which are just the entries of
the vector $b(Z)$) are the output of an MLP which takes in $X$. We learn the
weights of this neural net once up front from observed data $\mathcal{D}$ by
maximizing the log-likelihood of $\hat{p}(X\,|\,Z)$. In practice, this MLP has
$3$ hidden layers of sizes $(64,32,16)$ (first to last) and is trained with a
batch size of $512$ for $200$ epochs, with a learning rate of $0.01$. We use
ReLU activations.

Conditional normalizing flow. . We use an invertible (conditional) normalizing
flow to model the distribution of $X\,|\,Z$. Flows are a natural candidate for
modeling distributions, and in this case follow both the properties we want
from the modeling of $X\,|\,Z$. (i) Given $X,Z$, we can invert the
transformation to get $N$. (ii) We can sample from $X\,|\,Z=z$. We use the
Python library pyro Bingham et al.,  (2019).

All the experiments shown in this work use conditional normalizing flows due
to better performance and stable optimization. However, the _Gaussian
modeling_ is described to point out one more way in which prior assumptions on
functional dependencies could be incorporated in our generic framework. Once
again, the construction is identical for the leaky mediator setting. The only
difference is that there we model $M\,|\,X$ instead of $X\,|\,Z$.

###  E.3 Parameterizing $\mathcal{S}$

We consider $\theta$ to be $K$-dimensional and use the neural basis function
in all experiments unless specified otherwise. Equation 5 describes how we
parameterize our model. Here $\theta\,|\,N$ is parameterized by two MLPs with
parameters with weights denoted by $\eta_{0}$ and $\eta_{1}$. The weights of
this MLP constitute our main optimization parameters $\eta$ and we denote the
outputs by $\mu_{\eta_{0}}$ and $\Sigma_{\eta_{1}}$. $\mu_{\eta_{0}}$ has 2
hidden layers with size $(16,16)$ and $\Sigma_{\eta_{0}}$ has 2 hidden layers
with size $(32,32)$. To evaluate the objective function, we perform Monte
Carlo estimation $\operatorname{\mathbb{E}}_{N}[\mu_{\eta_{0}}(N)]$ with 500
samples for $N$ as in Equation 9.

The difference in the leaky mediator setting is that the parameterization of
$\theta$ has 2 noise inputs instead of 1, as described in Appendix B. The
remaining implementation choices remain unchanged.

Identifiable settings

Scalar treatment

(a) IV-lin-1d-weak-add. An identifiable scalar treatment with weak
confounding. We get tight bounds despite using the neural basis functions.

Multi-dimensional treatment

(b) IV-quad-2d-strong-add. Multi-dimension treatment using the neural basis
functions.

Figure 6: We see that our method is able to find tight bounds in identifiable
settings in the data-dense regions.

Strong instrument VS Weak instrument (Scalar treatment)

(a) iv-quad-1d-weak. This is a partially identifiable scalar treatment
setting. The bounds are particularly tight in the data-dense regions due to
the strong instrument.

(b) iv-quad-1d-strong. This is a partially identifiable scalar treatment
setting. The bounds are looser because the instrument is weak.

Figure 7: As expected, the bounds are tighter in the case of strong
instrument.

GAN comparison

(a) IV-quad-2d-weak.

(b) IV-quad-3d-weak.

Figure 8: We do note that in this instance, the bounds given by our method are
not as smooth in the data poor regions. However, our bounds are always valid,
while the GAN framework gives some invalid bounds in both cases above.

###  E.4 Response function choice

Our choice of response functions is characterized by the choice of $K$ basis
functions, as we saw in Equation 4. We primarily use the neural basis
functions described in Section 4.3. Note that a linear combination of basis
functions allows us to be arbitrarily expressive with out choice of family of
response functions, and we can choose In particular, we consider the following
options:

  1. 1.

_Neural basis functions:_ The basis function $\psi_{k}$ is the activation of
the $k$th neuron in the last hidden layer of an MLP which has been trained on
the observed data $\\{x_{i},y_{i}\\}_{i\in[n]}$ to learn $y$ given $x$.

In practice, we train a 3-hidden layer MLP with 64 neurons in the first two
layers, rectified linear units (ReLU) as activation functions and an MSE (mean
squared error) loss for $100$ epochs and a batch size of $512$ using Adam with
a learning rate of $0.01$. The size $K$ of the last hidden layer is equal to
the number of basis functions we wish to have in our family of response
functions.

An implicit assumption we have here is that the MLP can find a good
approximation of $y$ given $x$. This assumption is well supported by
theoretical and empirical results showing the approximation power of neural
networks (Shen et al., , 2021; Daubechies et al., , 2021), making it the
perfect candidate for multi-dimensional treatments.

  2. 2.

_Polynomials:_ For a multivariate input $x=(x_{1},x_{2}\ldots x_{d})$, a
response function can be considered to be a multivariate polynomial polynomial
in $x$. However, in this case we would have $2^{d}$ different basis functions
of degree up to $d$. This leads to the question of what basis functions to use
here so as to have enough expressive power, but not blow up the number of
basis functions and hence the dimension of the optimization parameter $\theta$
(see Equation 4). Due to this blowup in the choice of basis functions with the
dimension of the treatment in the case of polynomials, mainly use the neural
basis functions.

In the cases when we do use the polynomial basis, we restrict ourselves to
having upto quadratic terms in the basis function. For a $d$-dimensional
treatment, this amounts to having a $d(d+3)/2+1$ dimensional $\theta$.

  3. 3.

_Gaussian process basis functions (GP):_ We do not use this family of response
functions, but mention them as a possible option. There is some previous work
which considers the GP basis to define a family of functions (Kilbertus et
al., , 2020). In this approach a Gaussian process is fit to $K$ different sub-
samples $\\{(x_{i},y_{i})\\}_{i\in N^{\prime}}$ with $N^{\prime}\leq N$. A
single function is then sampled from each Gaussian process as the basis
functions $\psi_{k}$ for $k\in[K]$. Here, ‘sampling a function’ means to get
the evaluation of the function at several points in the treatment space and
then interpolate. While this is a reasonable approach when the treatment is
scalar, for a multidimensional treatment, this requires the interpolating a
multivariate function. The complexity and computational cost of such an
interpolation increases with the increase in dimensionality, while the
reliability of the interpolated function decreases at the same time. Also,
higher dimension require exponentially higher number of point evaluations to
get good interpolations. Due to these reasons, the GP basis approach is not
found suitable for higher dimensions. However, it can be a viable options for
scalar treatments.

###  E.5 Solving the optimization

The various MLPs described in the previous sections eventually lead to the
formulation of a constrained optimization problem that we have seen in
Equation 9.

We use the augmented Lagrangian method for inequality constraints to solve the
constrained optimization problem in Equation 9. The formulation is taken from
Section 17.3 in (Nocedal and Wright, , 2006).

We have seen in Section 3.5 that we have a total of $N\cdot L$ constraints. We
can think of $\hat{\phi}_{l}(x_{i},y_{i})$ as target values, estimated once up
front from observed data. We denote $\hat{\phi}_{l}(x_{i},y_{i})=B_{l,i}$. The
right-hand side $A_{l,i}$ is a function of the optimization parameter $\eta$
as seen in Section 3.4. For ease of notation, we “flatten” the indices $n$ and
$l$ into a single index $l\in[N\cdot L]$. We set the constraint slack to be
the same value for each constraint, so $\epsilon_{L}=\epsilon$.

Then our set of constraints is

$c_{l}(\eta):=\epsilon-|B_{l}-A_{l}(\eta)|\geq 0$

With this, the Lagrangian we aim to minimize with respect to $\eta$ can be
formulated as:

$\mathcal{L}(\eta,\lambda,\tau):=\pm o_{x^{\star}}(\eta)+\sum_{l=1}^{N\cdot
L}\xi(c_{l}(\eta),\lambda_{l},\tau)$ (14)

with

$\xi(c_{l}(\eta),\lambda_{l},\tau):=\begin{cases}-\lambda_{l}c_{l}(\eta)+\frac{\tau
c_{l}(\eta)^{2}}{2}&\text{if }\tau c_{l}(\eta)\leq\lambda_{l},\\\
-\frac{\lambda_{l}^{2}}{2\tau}&\text{otherwise},\end{cases}$

where $-$/$+$ is used for the upper/lower bound. $\tau$ is increased
throughout the optimization procedure and is seen as a temperature parameter.

Given an approximate minimum $\eta$ of this subproblem, we then update
$\lambda$ and $\tau$ according to
$\lambda_{l}\leftarrow\max\\{0,\lambda_{l}-\tau c_{l}(\eta)\\}$ and
$\tau\leftarrow\alpha\cdot\tau$ for all $l\in[N\cdot L]$ and a fixed
$\alpha>1$. The overall strategy is to iterate between minimizing Equation 14
and updating $\lambda_{l}$ and $\tau$. We find empirical justification in our
experiments, where the approach reliably converges on a range of different
datasets. We summarize our proposed procedure in Algorithm 1.

$l_{2}$ norm VS $l_{\infty}$ norm

Figure 9: IV-quad-weak. Partially identifiable multi-dimension treatment. The
$l_{2}$ norm includes all the constraints without sampling.

We have already seen how the basis function is chosen, how we can write the
$A$ of the constraints in closed form, how the $B$ of the constraints is fixed
upfront using $\phi_{1}$ and $\phi_{2}$, and how the objective is estimated.
It now remains, to solve this optimization problem, which we solve using the
augmented lagrangian method as described in S. In practice we fix the initial
value of $\tau$ to be $\tau_{init}=10$, and the multiplier $\alpha$ to be $5$,
meaning that $\tau$ is updated by multiplying it by $1.08$ at each step. The
max value of $\tau$ is fixed at $\tau_{max}=10,000$. We use $30$ optimization
steps to find the approximate optimal $x_{k}$ at the $k^{th}$ round of the
optimization, and perform $150$ round of optimization (number of times
$\lambda$ is updated) for each value of $x^{\star}$ and each bound (upper and
lower). The optimization was performed using the Adam optimizer with a
learning rate of $0.001$. All of this again is through the use of auto-
differentiation in PyTorch Paszke et al.,  (2019), implemented in python.

###  E.6 Hyperparameter search

We did not perform an automated hyperparameter search for $\tau_{init}$,
$\tau_{max}$ and $\tau_{factor}=\alpha$ since we observed that the values
$\tau_{init}=10$, $\tau_{max}=100000$ and $\tau_{factor}=5$ worked reasonably
well in all our settings. We arrived at this value through trial and error. An
automated hyperparameter search can also be expected to improve the bounds.

###  E.7 Computational resources

We used a computing cluster provided by the university, for ease of
parallelization of experiments. We used 1 CPU, 8000MB of RAM for all the
$x^{\star}$ bounds for a single random seed for the method. To clarify, the
use of the cluster was only for the purpose of running the optimization
algorithm parallely for various different random seeds. A single run of the
algorithm (for getting bounds on multiple values of $x^{\star}$) runs
comfortably on a local machine with 16GB of RAM. No GPUs were used for the
experiments.

##  Appendix F Further Experiments

Here we show the performance of our method in a variety of settings.

###  F.1 Identifiable settings

21 shows the performance in the case of identifiable settings. We see that our
method is able to find tight bounds in identifiable settings in the data dense
regions. However, it is natural that the more expressive the response
function, the wider the bounds. This can potentially be used as a test of
identifiability as described in Section 5. Such a test can be seen as a rather
continuous measure of identifiability, giving an indication of being somewhere
between identifiable and partially identifiable based on the tightness of
bounds, rather than being able to distinct strictly between identifiable and
non-identifiable settings.

###  F.2 Strong instrument VS Weak instrument

Figure 7 compares the bounds in the cases of having a strong or weak
instrument, or equivalently having weak or strong confounding relative to the
instrument. As expected the bounds are looser when the confounding is
stronger.

###  F.3 Comparing constraint norms

In Section 3.4 we have formulated our constraints in terms of norms. All the
experiments so far have been using the $\|\cdot\|_{\infty,\infty}$ norm. Here
we show experiments using the $\|\cdot\|_{2,2}$ norm as a constraint. This
optimization procedure remains the same, but is potentially easier to solve
since there is only $1$ constraint now, but on the other hand this constraint
is looser than matching the moment of each data point (or a chosen subsample).
The results can be seen in 24.

We point out that our formulation does not restrict us to norm based
constraints. In fact given the generative model we have define in Section 3,
we could also use any measure of distance between distributions as a
constraint. However this would again lead to solving a bi-level optimization,
making the optimization harder.

###  F.4 More comparisons to GAN bounds

We show some further comparisons of our method to the GAN bounds here. We see
a similar trend where our bounds are generally smoother. The plots can be
found in Figure 8, Figure 5(a) and Figure 5(b).

##  Appendix G Limitations of Kilbertus et al.,  (2020)

Kilbertus et al.,  (2020) limit the number of constraints by estimating
$\operatorname{\mathbb{E}}[\phi_{l}(Y)\,|\,z_{i}]$ only for a small pre-
determined set of grid points $z_{i}$ around which they bin the observed
datapoints for empirical estimates. Thereby, they only constrain $p(Y\,|\,Z)$
instead of $p(Y\,|\,X,Z)$ and fundamentally limit themselves to low-
dimensional settings, since otherwise the number of grid points grows
exponentially with the dimension of $\mathcal{Z}$. In particular, their
proposed interpolation of empirical cumulative distribution functions to fix
$p(x\,|\,z)$ in their copula model that parameterizes the distribution over
$\theta$ virtually only allows one-dimensional instruments and treatments.
Note that the binning procedure for $Z$ also becomes problematic in the small
data regime, where sufficiently many points are needed in each bin to keep the
variance in the empirical mean and variance estimates low. As a consequence,
the number of gridpoints in $z$ must be carefully tuned depending on the
dataset size.

Instead, our stochastic subsampling approach works with a fixed number of
constraints for each update step independent of the size of the dataset and
the dimensionality of $Z$ and $X$. In addition, we encode the structural
assumptions into a graphical model for $X,Z,\theta$ that still allows for
flexible conditional density estimation techniques such as invertible flows
for individual components. The copula model had to rely on interpolated
empirical cumulative density function estimated from potentially few
datapoints within each bin. The Gaussian copula is also less flexible than our
proposal. Finally, Kilbertus et al.,  (2020) used Monte Carlo estimates both
for the objective as well as for each individual constraint, often using
different sample sizes for these estimates. Instead, we exploit the form of
response functions as linear combinations of basis functions and compute the
constraints in closed form, removing the additional variance from the
stochastic optimization procedure. This leads to fewer tunable parameters and
more robust convergence of the optimization.

##  Appendix H Ethical and social implications

The ethical implications of causal inference are numerous, especially when
machine learning models are used to make high-stakes decisions in areas like
finance, criminal justice, and healthcare. Several authors in recent years
have developed causal notions of algorithmic fairness to acknowledge the
structural links between protected attributes and socially significant
outcomes Kusner et al.,  (2017); Kilbertus et al.,  (2017); Wu et al., 2019b ;
Kilbertus et al.,  (2019). More reliable methods for bounding causal effects
could help quantify the counterfactual fairness of individual decisions when
identifiability is impossible. By the same token, these bounds may
misrepresent treatment effects when untestable assumptions are not satisfied.
In the worst case, an agent may design a DAG in an adversarial manner to
ensure some desired result. This can be mitigated by public disclosure of all
structural assumptions, so that regulators and data subjects may critically
evaluate decision procedures.