# Causal Inference using Gaussian Processes
with Structured Latent Confounders

Sam Witty  Kenta Takatsu  David Jensen  Vikash Mansinghka

###### Abstract

Latent confounders—unobserved variables that influence both treatment and
outcome—can bias estimates of causal effects. In some cases, these confounders
are shared across observations, e.g. all students taking a course are
influenced by the course’s difficulty in addition to any educational
interventions they receive individually. This paper shows how to
semiparametrically model latent confounders that have this structure and
thereby improve estimates of causal effects. The key innovations are a
hierarchical Bayesian model, Gaussian processes with structured latent
confounders (GP-SLC), and a Monte Carlo inference algorithm for this model
based on elliptical slice sampling. GP-SLC provides principled Bayesian
uncertainty estimates of individual treatment effect with minimal assumptions
about the functional forms relating confounders, covariates, treatment, and
outcome. Finally, this paper shows GP-SLC is competitive with or more accurate
than widely used causal inference techniques on three benchmark datasets,
including the Infant Health and Development Program and a dataset showing the
effect of changing temperatures on state-wide energy consumption across New
England.

Machine Learning, ICML

##  1 Introduction

Multiple causal models can be observationally equivalent, i.e., they induce
the same likelihoods for observed data, while producing different estimates of
the effects of an intervention. This observational equivalance between causal
models is the basis for the colloquial expression “correlation does not imply
causation.” Distinguishing between causal models, and estimating the effects
of interventions, requires untestable assumptions about causal structure.

One such common assumption is unconfoundedness (Imbens & Rubin, 2015), i.e.,
that there exist no latent variables that influence both treatment and
outcome. This assumption enables the unique identification of interventional
distributions from the joint distribution over observed variables (Pearl,
2009) and reduces causal inference to probabilistic estimation. Unfortunately,
assuming unconfoundedness is often unreasonable in real observational settings
(Shadish et al., 2008). However, it may be more reasonable to assume
uncounfoundedness for a subset of data instances that are known to share a
common structure.

Symbol | Description | Entity
---|---|---
$\textbf{u}_{o}$ | Confounders | Object
$\mathrm{\textbf{x}}_{i}$ | Covariates | Instance
$\textbf{t}_{i}$ | Treatment | Instance
$\textbf{y}_{i}$ | Outcome | Instance
(a) Variable descriptions.

(b) Causal graph for GP-SLC.

$\displaystyle f_{u}$ $\displaystyle\sim GP(0,k_{u})$ $\displaystyle f_{x}$
$\displaystyle\sim GP(0,k_{x})$ $\displaystyle f_{t}$ $\displaystyle\sim
GP(0,k_{t})$ $\displaystyle f_{y}$ $\displaystyle\sim GP(0,k_{y})$

$\displaystyle\textbf{u}_{o=1...N_{O}}$
$\displaystyle=f_{u}(\boldsymbol{\epsilon}_{u_{o}})$
$\displaystyle\mathrm{\textbf{x}}_{i=1...N_{I}}$
$\displaystyle=f_{x}(\textbf{u}_{o=Pa(i)},\boldsymbol{\epsilon}_{x_{i}})$
$\displaystyle\textbf{t}_{i=1...N_{I}}$
$\displaystyle=f_{t}(\textbf{u}_{o=Pa(i)},\mathrm{\textbf{x}}_{i},\boldsymbol{\epsilon}_{t_{i}})$
$\displaystyle\textbf{y}_{i=1...N_{I}}$
$\displaystyle=f_{y}(\textbf{u}_{o=Pa(i)},\mathrm{\textbf{x}}_{i},\textbf{t}_{i},\boldsymbol{\epsilon}_{t_{i}})$

(c) Prior and causal functions for GP-SLC.

(d) Example grounding of the structural causal model in (b) and (c). Latent
confounders are shared within objects.

(e) Treatment, covariates, and inferred object-level confounders for instances
in (d). Color = $o$. Size = x.

(f) Kernel covariance matrix over observed ($\textbf{y}_{i}$) and
counterfactual ($\textbf{y}_{1,t_{*}}$) outcomes for instances in (e). Dark
$>$ light.

Figure 1: Model summary. GP-SLC (a-c) is a Gaussian process model for causal
inference in settings where object-level latent confounders, u, influence
instance-level observed covariates, x, treatment, t, and outcome, y, random
variables. For a given grounding (d), the outcome kernel function, $k_{y}$,
applied to treatment, covariates, and inferred confounders (e) induces the
covariance between observed and counterfactual outcomes (f). Instances
belonging to the same object always have the same inferred latent u. In this
example, the counterfactual outcome $\textbf{y}_{1,t_{*}}$ has high covariance
with factual outcomes $\textbf{y}_{1}$ and $\textbf{y}_{2}$.
$\textbf{y}_{1,t_{*}}$ has low, but non-zero, covariance with $\textbf{y}_{4}$
because $\textbf{u}_{Pa(1)}\not\approx\textbf{u}_{Pa(4)}$, despite the fact
that $\textbf{t}_{*}\approx\textbf{t}_{4}$ and
$\mathrm{\textbf{x}}_{1}\approx\mathrm{\textbf{x}}_{4}$.

For example, suppose an educator proposes a new policy of holding back poor
performing kindergarten students (Hong & Raudenbush, 2006; Hong & Yu, 2008)
with the intention of increasing their future academic performance. To
estimate the effect of this policy change, they gather data on student
retention and education outcomes from a national database. Here, the
unconfoundedness assumption is not justified, as the schools’ retention
policies are likely to be influenced by local economic conditions, which may
also influence student outcomes through other causal mechanisms, such as the
availability of educational resources. However, the assumption may be
justified when considering only students within a particular school, as this
subset of students are similarly influenced by local economic conditions. In
other words, statistical relationships within a school are less likely to be
biased by latent confounders than are statistical relationships across the
entire population.

In this paper, we present Gaussian processes with structured latent
confounders (GP-SLC), a novel Bayesian nonparametric approach to causal
inference with hierarchical data. The key innovation behind GP-SLC is to place
Gaussian process priors over functions in a hierarchical structural causal
model, bringing the flexibility of Gaussian process models to a wide variety
of practical causal inference techniques. GP-SLC naturally handles binary and
continuous treatments and requires minimal assumptions about functional
relationships between latent confounders, observed covariates, treatment, and
outcomes. See Figure 1 for an overview on how GP-SLC estimates counterfactual
outcomes from data.

##  2 Background

###  2.1 Object Conditioning

Recent work has studied how the analytical procedure of partitioning data
based on a known object hierarchy (e.g. students belonging to the same school)
relates to the syntax and semantics of causal graphical models (Jensen et al.,
2019). This work concludes that conditioning on the identify of objects
(referred to as object conditioning) is distinct from existing notions of
conditioning on the values of variables. Importantly, object conditioning
constrains a set of latent variables to be identical across a set of
instances, but does not constrain the particular value of those variables.
Furthermore, the statistical implications of object conditioning differ from
those of variable conditioning in that object conditioning does not induce
collider bias when variables on the object are caused jointly by treatment and
outcome.

Partitioning hierarchical data in this way is the key analytical procedure for
a variety of practical causal inference techniques, including within-subjects
designs (Loftus & Masson, 1994), difference-in-differences designs (Shadish et
al., 2008), longitudinal studies (Liang & Zeger, 1986), twin studies (Boomsma
et al., 2002), and multi-level-modeling (Gelman & Hill, 2006). As in the
student retention example, these techniques take advantage of background
knowledge about which instances (students) belong to which objects (schools)
to mitigate the biasing effects of latent confounders. However, these methods
typically rely on simple parametric assumptions, such as linear functional
dependencies. These parametric assumptions are often unjustified in real
domains, leading to poor estimates of causal effect.

We employ the idea of object conditioning directly in the GP-SLC model,
constraining the joint distribution over individuals’ latent confounders
instead of treating object identity as a covariate in and of itself. By
explicitly performing inference over object-level latent confounders, GP-SLC’s
estimates of counterfactual outcomes in one object are informed by observed
outcomes in another. Sharing information between objects in this way is
particularly valuable when each object contains few observed instances, as we
show in Section 6.

###  2.2 Causal Inference with Latent Confounders

Latent confounders—unobserved variables that cause both treatment and
outcome—bias estimates of treatment effect. However, the effect can be
adjusted for with additional background knowledge, such as that a latent
confounder influences an observed proxy variable (Kuroki & Pearl, 2014; Miao
et al., 2018). Similarly, recent work indicates that latent confounders can be
adjusted for if they cause multiple candidate treatment variables (Wang &
Blei, 2019).

GP-SLC is similar to these approaches, in that it leverages additional
background knowledge to adjust for latent confounders. However, unlike prior
work using generative models for causal inference with latent confounders, it
leverages known hierarchical structure to identify causal effects.

###  2.3 Gaussian Process Models

Gaussian process models are a flexible technique for probabilistic modeling.
Specifically, a Gaussian process is a distribution over deterministic
functions, $\textbf{y}=f(\mathrm{\textbf{x}})$, $f\sim GP(m,k)$, which is
fully specified by its mean function, $m(\mathrm{\textbf{x}})$ and covariance
function, $k(\mathrm{\textbf{x}},\mathrm{\textbf{x}}^{\prime})$, which we will
refer to as the kernel function (Rasmussen, 2003). By definition, any finite
collection of draws from a Gaussian process prior are jointly Gaussian
distributed, $Y\sim\mathcal{N}(\mu,\Sigma)$, where
$\mu_{i}=m(\mathrm{\textbf{x}}_{i})$ and
$\Sigma_{i,{i^{\prime}}}=k(\mathrm{\textbf{x}}_{i},\mathrm{\textbf{x}}_{i^{\prime}})$.
We denote such covariance matrices as $K(X,X)$, where $X$ and $Y$ are matrices
of all $\mathrm{\textbf{x}}_{i}$ and $\textbf{y}_{i}$ respectively. It is
common to set the prior mean function to $m(\mathrm{\textbf{x}})=0$, which we
do in GP-SLC.

This identity is useful for two reasons: (i) it provides an explicit
likelihood, which can be used to perform inference over latent variables
(Lawrence, 2004; Titsias & Lawrence, 2010); and (ii) it enables closed-form
out-of-sample probabilistic prediction (Rasmussen, 2003). We take advantage of
both of these characteristics in GP-SLC, performing approximate inference over
latent confounders in Section 4.1 and predicting counterfactual outcomes in
Section 4.2.

###  2.4 Structural Causal Models

GP-SLC can be thought of as a Bayesian nonparametric prior over functions in a
structural causal model (SCM). SCMs provide a syntax and semantics for
reasoning about interventional and counterfactual distributions in a system of
random variables (Pearl, 2009). Counterfactuals (Pearl, 2011)—answers to what-
if questions—are expressed in probability notation as
$P(\textbf{y}_{\textbf{t}_{*}}|\textbf{y},\textbf{t})$, where
$\textbf{y}_{\textbf{t}_{*}}$ is the counterfactual outcome under intervention
$do(\textbf{t}=\textbf{t}_{*})$, y is the observed outcome, and t is the
observed treatment. In our education example, we may be interested in the
counterfactual, “given that the student was not retained in kindergarten and
they performed poorly in high school, how would they have performed if they
had been retained?” Answering these kinds of counterfactual queries involves:
(i) computing the conditional distribution over latent variables and exogenous
noise given observed evidence; (ii) applying the intervention to the
structural causal model; and (iii) recomputing the distribution over the
outcome variable(s) using the modified structural causal model. We apply this
procedure to estimate counterfactual outcomes using GP-SLC in Section 3.

##  3 Gaussian Processes with Structured Latent Confounders

Consider the common scenario where there are $N_{O}$ object-level latent
confounders $(\mathit{\textbf{u}_{o}\in\mathbb{R}^{N_{U}}})$ that influence
$N_{I}$ instances of observed treatment $(\textbf{t}_{i}\in\mathbb{R})$,
covariates $(\mathrm{\textbf{x}}_{i}\in\mathbb{R}^{N_{X}})$, and outcomes
($\textbf{y}_{i}\in\mathbb{R}$). We can describe this scenario as a structural
causal model, where the particular functions relating u, x, t, and y are given
by the following:

$\displaystyle\begin{aligned}
\textbf{u}_{o=1...N_{O}}&=f_{u}(\boldsymbol{\epsilon}_{u_{o}})\\\
\mathrm{\textbf{x}}_{i=1...N_{I}}&=f_{x}(\textbf{u}_{o=Pa(i)},\boldsymbol{\epsilon}_{x_{i}})\\\
\textbf{t}_{i=1...N_{I}}&=f_{t}(\textbf{u}_{o=Pa(i)},\mathrm{\textbf{x}}_{i},\boldsymbol{\epsilon}_{t_{i}})\\\
\textbf{y}_{i=1...N_{I}}&=f_{y}(\textbf{u}_{o=Pa(i)},\mathrm{\textbf{x}}_{i},\textbf{t}_{i},\boldsymbol{\epsilon}_{y_{i}}).\end{aligned}$
(1)

If all instances belong to the same object ($N_{O}=1$) the structural causal
model in Equation 1 reduces to the standard propositional case and the latent
u will not bias estimated counterfactual outcomes. However, if we wish to
estimate counterfactual outcomes using instances from multiple objects
($N_{I}>N_{O}>1$), u’s influence on t and y would appear to render
counterfactual queries unidentifiable (Pearl, 2009). To address this problem
of identifiability, GP-SLC places a mean-zero Gaussian process prior over each
function in the structural causal model in Equation 1, with kernel functions
$k_{x}$, $k_{t}$, and $k_{y}$ respectively as follows:

$\displaystyle\begin{aligned} f_{x}&\sim GP(0,k_{x})&f_{t}&\sim
GP(0,k_{t})&f_{y}&\sim GP(0,k_{y}).\end{aligned}$

The particular choice of each kernel function plays an important role in the
prior over functions, and by extension the conditional distribution over
counterfactual outcomes. We use a radial basis function (RBF) kernel with
automatic relevance determination (ARD) (Neal, 2012) and additive Gaussian
exogenous noise for each Gaussian process prior. Each kernel is parameterized
by a set of kernel lengthscales, $\lambda$, scaling factors, $\sigma^{2}$, and
exogenous noise variances $\sigma^{2}_{\epsilon}$. We assume $f_{u}$ is the
identity function. We refer to the noise-free component of each kernel
function as $k^{\prime}$, e.g.
$k_{t}=k^{\prime}_{t}([\textbf{u}_{o=Pa(i)},\mathrm{\textbf{x}}_{i}],[\textbf{u}_{o^{\prime}=Pa(i^{\prime})},\mathrm{\textbf{x}}_{i^{\prime}}])+\sigma^{2}_{\epsilon_{y}}\delta_{i,i^{\prime}}$,
where $\sigma^{2}_{\epsilon_{y}}$ is the exogenous noise variance,
$\delta_{i,i^{\prime}}$ is the Dirac-delta function at $i^{\prime}=i$, and
$k^{\prime}_{t}$ is the ARD kernel. See the supplementary materials for
detailed descriptions of these kernels.

In addition to placing Gaussian process priors on the functions in the
structural causal model in Equation 1, we also place inverse-gamma priors,
$P(\theta)=\gamma\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}(\theta;\alpha_{\theta},\beta_{\theta})$
on each $\theta\in\Theta$, where $\Theta$ is the set of all kernel
lengthscales, scaling factors, and exogenous noise variances. In Section 4.1
we show how to perform approximate posterior inference on $\Theta$.

###  3.1 Conditional Density

As $f_{y},f_{t},\text{ and }f_{x}$ are all drawn from Gaussian process priors,
$P(Y|T,X,U,\Theta),P(T|X,U,\Theta),\text{ and }P(X|U,\Theta)$ are all
multivariate Gaussian distributions with mean zero and covariance given by
their respective kernel covariance matrices. For example, $P(T|X,U,\Theta)=$
$\mathcal{N}(T;0,K_{t}([U,X],[U,X]))$, where
$K_{t}([U,X],[U,X])_{i,i^{\prime}}$
$=k_{t}([\textbf{u}_{o=Pa(i)},\mathrm{\textbf{x}}_{i}],[\textbf{u}_{o^{\prime}=Pa(i^{\prime})},\mathrm{\textbf{x}}_{i^{\prime}}])$.
As $\textbf{u}_{o}$ is given by the identity function of exogenous Gaussian
noise, $P(\textbf{u}_{o}|\Theta)$
$=\mathcal{N}(\textbf{u}_{o};0,\sigma^{2}_{\epsilon_{u}}\textit{I}_{N_{U}})$.
Therefore, the joint density is given by the following, which we use in
Algorithms 2 and 3:

Algorithm 1 Individual Treatment Effect Estimation

1: Input:

2: Intervention assignment: $\textbf{t}_{*}$

3: Data: $Y,T,X$

4: Prior hyperparameters: $\alpha_{\theta\in\Theta},\beta_{\theta\in\Theta}$

5: Inference parameters:
$\textrm{N}_{\textrm{Outer}},\textrm{N}_{\textrm{MH}},\textrm{N}_{\textrm{ES}},\textrm{drift}_{\theta\in\Theta}$

6: Procedure:

7:
$\theta\sim\gamma\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}(\alpha_{\theta},\beta_{\theta}),\forall\theta\in\Theta$
$\triangleright$ Prior sample

8: $\textbf{u}_{o}\sim\mathcal{N}(0,\sigma^{2}_{U}\textit{I}_{N_{U}}),\forall
o=1...N_{O}$ $\triangleright$ Prior sample

9: $\textit{ITESamples}\leftarrow\\{\hbox{}\\}$

10: for $l=1$ to $n_{Outer}$ do

11: $\Theta\leftarrow\text{HyperparameterUpdate}(...)$ $\triangleright$
Algorithm 2

12: $U\leftarrow\text{ConfounderUpdate}(...)$ $\triangleright$ Algorithm 3

13:
$\textbf{w}_{i}\leftarrow[\textbf{t}_{i},\mathrm{\textbf{x}}_{i},\textbf{u}_{o=pa(i)}],\forall
i\in 1...N_{I}$

14:
$\textbf{w}_{i,*}\leftarrow[\textbf{t}_{*},\mathrm{\textbf{x}}_{i},\textbf{u}_{o=pa(i)}],\forall
i\in 1...N_{I}$

15:
$\mu_{\textit{ITE}}\leftarrow(K^{\prime}(W,W_{*})\text{-}K^{\prime}(W,W))K(W,W)\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}Y$

16: $\textit{ITE}\sim\mathcal{N}(\mu_{\textit{ITE}},\Sigma_{\textit{ITE}})$
$\triangleright$ See Supplement for $\Sigma_{\textit{ITE}}$

17: $\textit{ITESamples}\leftarrow\textit{ITE}\cup\textit{ITESamples}$

18: return ITESamples

$\displaystyle\begin{aligned}
P(Y,T,X,U,\Theta){}=&P(Y|T,X,U,\Theta)P(X|U,\Theta)\\\
&P(T|X,U,\Theta)\prod_{o=1...N_{O}}P(\textbf{u}_{o}|\Theta)P(\Theta).\end{aligned}$

By placing Gaussian process priors over each function in the hierarchical
structural model, we encode our assumptions about which configurations of
observed and latent variables are reasonable a-priori. Using a radial basis
function kernel, we assume that if two objects have similar object-level
latent confounders, they are likely to induce similar distributions over
observed covariates, treatment, and outcome. Placing higher density on smooth
structural causal functions in this way enables inference over object-level
confounders.

##  4 Estimating Treatment Effects

In this section we describe how to estimate the individual treatment effect,
$\textit{ITE}_{i,t_{*}}=\textbf{y}_{i,t_{*}}-\textbf{y}_{i}$, the difference
between observed and counterfactual outcomes for the $i$th instance. Standard
aggregate measures of causal effect, such as the sample average treatment
effect, $\textit{SATE}_{t_{*}}=\frac{1}{N_{I}}\sum_{i}\textit{ITE}_{i,t_{*}}$,
can be derived from the individual treatment effect. We use
$\textit{ITE}_{t_{*}}$ to denote the vector of individual treatment effects
for the intervention $do(\textbf{t}_{i}=\textbf{t}_{*})$ applied uniformly to
each instance $i$, although the estimation procedure can be easily applied to
any arbitrary set of intervention assignments.

First, note that when exogenous noise is additive in $f_{y}$, i.e
$f_{y}(\textbf{u}_{o=Pa(i)},\mathrm{\textbf{x}}_{i},\textbf{t}_{i},\boldsymbol{\epsilon}_{t_{i}})=$
$f^{\prime}_{y}(\textbf{u}_{o=Pa(i)},\mathrm{\textbf{x}}_{i},\textbf{t}_{i})+\boldsymbol{\epsilon}_{t_{i}}$,
as in the GP-SLC model, individual treatment effect is given by the difference
between noise-free functions
$\textit{ITE}_{i,t_{*}}=f^{\prime}_{y}(\textbf{u}_{o=Pa(i)},\mathrm{\textbf{x}}_{i},\textbf{t}_{*})-f^{\prime}_{y}(\textbf{u}_{o=Pa(i)},\mathrm{\textbf{x}}_{i},\textbf{t}_{i})$.
We denote the outcome of these noise-free functions as
$\textbf{y}^{\prime}_{i,t_{*}}$ and $\textbf{y}^{\prime}_{i}$, and the vector
of outcomes as $Y^{\prime}_{t_{*}}$ and $Y^{\prime}$ respectively.111Noise-
free prediction is often denoted as $f$ in Gaussian process regression models.
We avoid this notation to avoid confusion with functions in the structural
causal model. As $U\cup X$ blocks all backdoor paths from $T$ to $Y$, we have
that the distribution over individual treatment effects is given by the
following expression (Pearl, 2009):

$\displaystyle\begin{aligned}
&P(\textit{ITE}_{t_{*}}|Y,T,X)=P(Y^{\prime}_{t_{*}}-Y^{\prime}|Y,T,X)\\\
&=\int
P(Y^{\prime}_{*}-Y^{\prime}|T_{*},Y,T,X,U,\Theta)P(U,\Theta|Y,T,X)dUd\Theta.\end{aligned}$

This equation directly informs our hybrid procedure for estimating
counterfactual outcomes shown in Algorithm 1, (i) generate approximate samples
from the posterior $\hat{U},\hat{\Theta}\sim P(U,\Theta|Y,T,X)$ and (ii) for
each posterior sample compute the conditional distribution
$(Y^{\prime}_{*}-Y^{\prime}|T_{*},Y,T,X,\hat{U},\hat{\Theta})$ in closed-form,
taking advantage of Gaussian closure under conditioning and subtraction. As
the posterior distribution $P(U,\Theta|Y,T,X)$ is intractable for non-trival
kernels, we turn to Monte Carlo approximate inference techniques.

###  4.1 Approximate Inference: $U$ and $\Theta$

Because we assume that our structural functions were drawn from Gaussian
Processes, which provide a closed-form expression for the conditional density
of the data, we are able to use standard likelihood-based approximate
inference techniques. In our experiments, we approximate this posterior
distribution using elliptical slice sampling (Murray et al., 2010) for the
latent confounder, $U$, and random walk Metropolis Hastings (Hastings, 1970)
on all kernel hyperparameters and exogenous noise variances, $\Theta$. Psuedo-
code implementations are presented in Algorithms 2 and 3.

###  4.2 Exact Inference: $Y^{\prime}_{*}-Y^{\prime}$

To estimate $P(Y^{\prime}_{*}-Y^{\prime}|T_{*},Y,T,X,U,\Theta)$, we extend the
Gaussian process model over in-sample and out-of-sample outcomes (Rasmussen,
2003). Using the shorthand
$\textbf{w}_{i}=[\textbf{t}_{i},\mathrm{\textbf{x}}_{i},\textbf{u}_{o=pa(i)}]$
and
$\textbf{w}_{i,*}=[\textbf{t}_{*},\mathrm{\textbf{x}}_{i},\textbf{u}_{o=pa(i)}]$,
the joint distribution over observed outcomes, $Y$, noise-free outcomes for
each observed instance, $Y^{\prime}$, and noise-free counterfactual outcomes,
$Y^{\prime}_{*}$ conditioned on observed treatments, $T$, covariates, $X$,
inferred confounders, $U$, and kernel hyperparameters, $\Theta$, is Gaussian
distributed as follows, where
$K(W,W)=K^{\prime}(W,W)+\sigma^{2}_{Y}\textit{I}_{N_{I}}$ and
$K^{\prime}(W,W)$ is the kernel matrix of $k^{\prime}_{y}$ given $\Theta$:

$\displaystyle\Bigg{(}\begin{bmatrix}[l]Y\\\ Y^{\prime}\\\
Y^{\prime}_{*}\end{bmatrix}$ $\displaystyle|T_{*},T,X,U,\Theta\Bigg{)}$
$\displaystyle\sim\mathcal{N}$
$\displaystyle\Bigg{(}0,\begin{bmatrix}[l]K(W,W)&K^{\prime}(W,W)&K^{\prime}(W,W_{*})\\\
K^{\prime}(W,W)&K^{\prime}(W,W)&K^{\prime}(W,W_{*})\\\
K^{\prime}(W_{*},W)&K^{\prime}(W_{*},W)&K^{\prime}(W_{*},W_{*})\end{bmatrix}\Bigg{)}.$

As Gaussian distributions are closed under conditioning and subtraction, we
have that $(Y^{\prime}_{*}-Y^{\prime}|T_{*},Y,T,X,U,\Theta)$ is also jointly
Gaussian distributed as follows, where
$\mu_{\textit{ITE}}=(K^{\prime}(W,W_{*})-K^{\prime}(W,W))K(W,W)\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}Y$:

$\displaystyle\begin{aligned}
(Y^{\prime}_{*}-Y^{\prime}|&T_{*},Y,T,X,U,\Theta)\sim\mathcal{N}(\mu_{\textit{ITE}},\Sigma_{\textit{ITE}})\end{aligned}$
(2)

See the supplementary materials for details and for a closed-form expression
for $\Sigma_{\textit{ITE}}$.

Algorithm 2 Hyperparameter Update - Random Walk MH

1: Input:

2: Data: $Y,T,X$

3: Posterior sample: $U,\Theta$

4: Prior hyperparameters: $\alpha_{\theta\in\Theta},\beta_{\theta\in\Theta}$

5: Inference parameters:
$\textrm{N}_{\textrm{MH}},\textrm{drift}_{\theta\in\Theta}$

6: Procedure:

7: for $j=1$ to $\textrm{N}_{\textrm{MH}}$ do

8: for $\theta\in\Theta$ do

9: $\alpha_{\theta^{\prime}}\leftarrow\theta^{2}/\textrm{drift}_{\theta}$

10: $\beta_{\theta^{\prime}}\leftarrow\theta(\alpha_{\theta^{\prime}}-1)$

11:
$\theta^{\prime}\sim\gamma\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}(\alpha_{\theta^{\prime}},\beta_{\theta^{\prime}})$

12: $\alpha_{\theta^{\prime\prime}}\leftarrow\theta^{\prime
2}/\textrm{drift}_{\theta^{\prime}}$

13:
$\beta_{\theta^{\prime\prime}}\leftarrow\theta^{\prime}(\alpha_{\theta^{\prime\prime}}-1)$

14: $\Theta^{\prime}\leftarrow\Theta\setminus\theta\cup\theta^{\prime}$

15:
$A\leftarrow\dfrac{P(Y,T,X,U,\Theta^{\prime})}{P(Y,T,X,U,\Theta)}\dfrac{\gamma\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}(\theta^{\prime};\alpha_{\theta^{\prime}},\beta_{\theta^{\prime}})}{\gamma\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}(\theta;\alpha_{\theta^{\prime\prime}},\beta_{\theta^{\prime\prime}})}$

16: $\eta\sim\text{Uniform}(0,1)$

17: if $\eta>\min(A,1)$ then

18: $\Theta=\Theta^{\prime}$

19: return $\Theta$

Algorithm 3 Confounder Update - Elliptical Slice Sampling

1: Input:

2: Data: $Y,T,X$

3: Posterior sample: $U,\Theta$

4: Inference parameter: $\textrm{N}_{\textrm{ES}}$

5: Procedure:

6: for $j=1$ to $\textrm{N}_{\textrm{ES}}$ do

7: for $k=1$ to $\textrm{N}_{\textrm{U}}$ do

8: $\textrm{done}\leftarrow\textrm{False}$

9: $\nu\sim\mathcal{N}(0,\sigma^{2}_{\epsilon_{U}}\textit{I})$

10: $y\sim\text{Uniform}(0,P(Y,T,X,U,\Theta))$

11: $\phi\sim\text{Uniform}(0,2\pi)$

12: $[\phi_{\textrm{min}},\phi_{\textrm{max}}]\leftarrow[\phi-2\pi,\phi]$

13: while not done do

14: $U^{\prime}_{*,k}\leftarrow U_{*,k}\cos{\phi}+\nu\sin{\phi}$

15: if $P(Y,T,X,U,\Theta)>y$ then

16: $\textbf{u}_{k}\leftarrow\textbf{u}^{\prime}_{i}$

17: $\textrm{done}\leftarrow\textrm{True}$

18: else

19: if$\;\phi<0$ then$\;\phi_{\textrm{min}}\leftarrow\phi$
else$\;\phi_{\textrm{max}}\leftarrow\phi$

20: $\phi\sim\text{Uniform}(\phi_{\textrm{min}},\phi_{\textrm{max}})$

21: return $U^{\prime}$

##  5 Asymptotic Posterior Consistency

In the special case where each RBF kernel in the GP-SLC model is replaced with
a linear kernel, $k(A,A^{\prime})=A\cdot A^{\prime}$, shared confounding among
instances enables asymptotically consistent estimates of individual treatment
effect. This is contrasted with the propositional setting (i.e. $N_{O}=N_{I}$)
which does not lead to asymptotically consistent counterfactual estimation.
Informally, a continuous random variable $\psi$ is asymptotically consistent
if its posterior $P(\psi|data)$ approaches a Dirac-delta distribution at some
point $\psi^{\prime}$, regardless of the prior $P(\psi)$. We present proofs of
Proposition 5.1 and Theorems 5.2 and 5.3 in the supplementary materials.

######  Proposition 5.1.

When $N_{O}=N_{I}$, $\textit{ITE}_{t_{*}}$ is not asymptotically consistent
$\forall t_{*}\in\mathbb{R}$.

######  Theorem 5.2.

Assume there exists an object $o$ that is the parent of $n$ instances,
$I^{\prime}=\\{i^{\prime}_{1},...,i^{\prime}_{n}\\}$. Then
$\textit{ITE}_{t_{*}}$ is asymptotically consistent as $n$ approaches
$\infty,\forall t_{*}\in\mathbb{R}$.

######  Theorem 5.3.

Assume there exist $n$ objects $\mathbb{O}=\\{o_{1},...,o_{n}\\}$, each of
which are the unique parents of $k\geq 2$ instances
$I^{\prime}_{o}=\\{i^{\prime}_{o,1},...,i^{\prime}_{o,k_{o}}\\}$. Then
$\textit{ITE}_{t_{*}}$ is asymptotically consistent as $n$ approaches
$\infty,\forall t_{*}\in\mathbb{R}$.

(a) Original data.

(b) Unbiased sampling.

(c) Biased sampling.

(d) Energy consumption (GWh)

(e) Mean squared error in estimated sample average treatment effect.

Figure 2: Process and results for New England energy consumption benchmark. We
sample hotter days with higher probability for states with higher daily energy
consumption (a-d). Sampling in this way simulates confounding, creating an
observational relationship (consumption is signicantly higher in hotter days)
that differs from the causal relationship (low or high temperature causes a
moderate increase in energy consumption). GP-SLC (this paper) produces
accurate estimates of counterfactual outcomes, despite this confounding bias
(e). For baselines that ignore hierarchical structure (GP-NoObj and GP-
NoConf), accuracy decreases significantly with increasing confounding bias.
Results are normalized by the $\sqrt{\textrm{MSE}}$ of the GP-SLC model with
$\textrm{bias}=9^{\circ}F$ and 25 samples per state.

##  6 Experiments

Unlike associational models, which can be evaluated using accuracy on held-out
test data, causal models produce predictions about unobserved counterfactual
distributions. As a result, effective evaluation of causal models requires
different methods (Gentzel et al., 2019). We evaluate the GP-SLC model using
three benchmarks with known counterfactual outcomes. In Section 6.1, we
evaluate GP-SLC using a fully synthetic hierarchical data generating process.
In Section 6.2 we modify the Infant Health and Development Program (IHDP)
benchmark (Hill, 2011) to include hierarchical structure and latent
confounders. In Section 6.3 we introduce and evaluate on a new benchmark task
for observational causal inference with hierarchical data, predicting the
effect of changes in temperature on state-wide electric energy consumption in
New England (NEEC).

We implement the GP-SLC model using Gen (Cusumano-Towner et al., 2019), a
probabilistic programming language with programmable inference. Except where
otherwise specified we set $N_{U}=3$ and $\alpha_{\theta}=\beta_{\theta}=4$
for each inverse gamma prior over kernel hyperparameters and exogenous noise
variance. We estimate individual treatment effects using Algorithm 1, with
$\textrm{N}_{\textrm{Outer}}=5000$, $\textrm{N}_{\textrm{MH}}=3$,
$\textrm{N}_{\textrm{ES}}=5$, and
$\textrm{drift}_{\theta}=0.5,\forall\theta\in\Theta$.

We compare the GP-SLC model against six baselines: a GP regression model that
ignores latent confounding variables (GP-NoConf), a GP-SLC model where each
instance is incorrectly assigned a single object (GP-NoObj), a seperate GP
regression model for each object (GP-PerObj), Bayesian additive regression
trees (BART) (Hill, 2011), a random slope and intercepts linear model (MLM 1),
and a random intercepts linear model (MLM 2) (Gelman, 2006). The Gaussian
process baselines are ablations of the full GP-SLC model, and use the same
kernels, priors over hyperparameters, and inference scheme. The BART baseline
uses the object identifier, $o$, as an additional covariate. See the
supplementary materials for additional details on baselines.

We use two evaluation metrics to evaluate GP-SLC and baselines, mean squared
error of the sample average treatment effect,
$\text{MSE}=\mathbb{E}_{t_{*}}[(\textit{SATE}^{*}_{t_{*}}-\textit{SATE}_{t_{*}})^{2}]$,
and precision in estimation of heterogenous effect (Hill, 2011),
$\text{PEHE}=\mathbb{E}_{t_{*}}[\sum_{i}^{N_{i}}(\textit{ITE}^{*}_{i,t_{*}}-\textit{ITE}_{i,t_{*}})^{2}/N_{i}]$,
where $\textit{ITE}^{*}_{i,t_{*}}$ and $\textit{SATE}^{*}_{T_{*}}$ are the
actual effects and $\textit{ITE}_{i,t_{*}}$ and $\textit{SATE}_{t_{*}}$ are
the predicted effects. For the synthetic benchmark, we average over 100
regular intervals between the 5th and 95th percentile of treatment assignment
in the observational data. For the NEEC benchmask, we average over
$\\{30,30.1,...,70^{\circ}{F}\\}$.

Figure 3: Comparison among methods on the New England energy consumption
benchmark. Above are GP-SLC and all baselines’ effect estimates on the NEEC
benchmark with $\text{bias}=9^{\circ}F$ and 25 samples per state. Green shaded
regions indicate $90\%$ credible intervals. GP-SLC effectively recovers the
effect of temperature on energy consumption, despite the latent confounding
introduced by biased sampling. The best performing baseline, GP-PerObj,
produces poor estimates of the effect of high temperatures in Rhode Island.

###  6.1 Synthetic Data

We evaluate GP-SLC and various baselines on two synthetic datasets with
hierarchically structured latent confounders, one with additive and one with
multiplicative treatment and outcome functions. Both synthetic datasets are
generated using three dimensional object-level confounders for 20 objects,
each of which contains 10 instances. Observed instance-level covariates are
generated as a linear function of object-level Gaussian distributed latent
confounders. Details for synthetic treatment and outcome functions are
presented in the supplementary materials, and evaluation results are shown in
Table 1. GP-SLC consistently matches and exceeds the counterfactual prediction
performance of the six baselines on synthetic data. Baselines that ignore
object structure (GP-NoConf and GP-NoObj) produce the least accurate
counterfactual predictions.

Model | Additive | Multiplicative
---|---|---
$\sqrt{\text{PEHE}}$ | $\sqrt{\text{MSE}}$ | $\sqrt{\text{PEHE}}$ | $\sqrt{\text{MSE}}$
GP-SLC | 1.0 | 1.0 | 1.0 | 1.0
GP-NoConf | 21.3 | 25.3 | 4.2 | 7.6
GP-NoObj | 22.2 | 27.0 | 4.5 | 8.1
GP-PerObj | 3.7 | 3.4 | 1.1 | 0.9
MLM1 | 1.2 | 1.02 | 2.4 | 2.9
MLM2 | 1.3 | 1.6 | 4.4 | 9.3
BART | 8.5 | 10.7 | 2.6 | 4.3
Table 1: Results on synthetic data with additive and multiplicative nonlinear
treatment and outcome functions. Scores are normalized by the score of GP-SLC.
Lower is better.

In addition to the synthetic experiments presented in Table 1, we tested the
behaviour of GP-SLC using two alternative synthetic data generating processes.
On the first, a linear structural data generating process with shared
confounding, GP-SLC produces comparable estimates to the multi-level model
baselines. On the second, in which each object shares a common effect of
treatment and outcome rather than a common cause, GP-SLC is not susceptible to
collider bias (Berkson, 1946; Elwert & Winship, 2014). This empirical finding
is consistent with recent theory on object conditioning (Jensen et al., 2019).

###  6.2 Infant Health and Development Program

The IHDP benchmark (Hill, 2011) uses real data for treatments (whether a child
recieves high-quality child care and home visits from a trained provider) and
covariates (birth weight, head circumference, etc.) from the 1992 Infant
Health and Development Program (Ramey et al., 1992) with a synthetic nonlinear
outcome function. We modify the IHDP benchmark to simulate hierarchically
structured data by randomly duplicating 30% of the data instances and
reassigning the duplicate’s treatment assignment to be the opposite of the
original instance. In order to introduce variation between duplicated
instances, we add noise to each individuals’ continuous covariates from a
$\mathcal{N}(0,\sigma^{2}_{j})$, where $\sigma^{2}_{j}$ is 5% of the $j$th
covariate’s marginal variance. We obscure the remaining 15 categorical
covariates, representing object-level latent confounding. Even though the 15
categorical covariates are obscured from the GP-SLC model, they are identical
across duplicates, unlike the observed covariates. We then generate observed
and counterfactual outcomes using the benchmark synthetic outcome function,
applied to treatment, modified covariates, and latent confounders. In this
setting, $Pa(i)=Pa(i^{\prime})$ if instance $i$ is a duplicate of instance
$i^{\prime}$ or vice versa. Although each duplicate’s treatment assignment is
deterministic, the overall relationship between treatment and outcome is still
confounded, as we only duplicate a subset of the original instances.

For the IHDP benchmark, which has binary treatment variables, we modify the
GP-SLC model by replacing the expression
$\textbf{t}_{i}=f_{t}(\textbf{u}_{o=Pa(i)},\mathrm{\textbf{x}}_{i},\boldsymbol{\epsilon}_{t_{i}})$
with the expressions
$\hat{\textbf{t}}_{i}=f_{\hat{\textbf{t}}}(\textbf{u}_{o=Pa(i}),X_{i},\boldsymbol{\epsilon}_{\hat{t}_{i}})$
and
$\textbf{t}_{i}\sim\text{Bernoulli}(\textit{expit}(\hat{\textbf{t}}_{i}))$. In
this setting, we use elliptical slice sampling to approximate the latent logit
probability of treatment, $\hat{\textbf{t}}$.

Given the small size of each object, we omit the GP-PerObj baseline model from
this evaluation. As the IHDP benchmark includes binary treatment variables we
compared against four additional baselines: balanced linear regression
(BalReg) and balanced neural nets (BALNN) (Johansson et al., 2016), targeted
maximum likelihood estimation with the superlearner (TMLE) (Van der Laan et
al., 2007), and inverse probability of treatment weighting with logistic
regression (IPTW) (Imbens & Rubin, 2015).

Results of the IHDP evaluation are presented in Table 2. GP-SLC matches and
exceeds the performance of other baselines when predicting the effect of
assigning treatment to individuals who were previously untreated. In this
setting, the linear models (MLM 1 and MLM 2) produce the least accurate
counterfactual predictions.

Model | Control | Treated
---|---|---
$\sqrt{\text{PEHE}}$ | $\sqrt{\text{MSE}}$ | $\sqrt{\text{PEHE}}$ | $\sqrt{\text{MSE}}$
GP-SLC | 1.0 | 1.0 | 1.0 | 1.0
GP-NoConf | 1.03 | 1.07 | 1.04 | 0.94
GP-NoObj | 1.11 | 1.02 | 0.82 | 1.08
MLM1 | 68.3 | 33.2 | 106.7 | 1028.4
MLM2 | 73.3 | 389.1 | 45.8 | 63.2
BART | 3.7 | 1.1 | 2.4 | 0.33
BALReg | 5.1 | 82.7 | 1.9 | 0.5
BALNN | 2.1 | 7.0 | 1.7 | 4.5
TMLE | n/a | 209.8 | n/a | 12.2
IPTW | n/a | 50.6 | n/a | 90.5
Table 2: Results on the modified infant health and development program
benchmark, shown seperately for treated and untreated individuals. Scores are
normalized by the score of GP-SLC. TMLE and IPTW do not estimate individual
treatment effects. Lower is better.

###  6.3 New England Energy Consumption

We introduce a new benchmark for estimating heterogenous effects in
hierarchically structured settings, predicting the effect of changing
temperature on state-wide electric energy consumption in New England. Unlike
the evaluation in Section 6.2, which includes real treatments, covariates, and
confounders and a synthetic outcome function, the New England energy
consumption (NEEC) benchmark preserves outcome functions from real quasi-
experimental data, and uses biased sampling to induce confounding.
Specifically, we generate data for the NEEC benchmark task using the New
England Independent Service Operator’s public records on hourly dry-bulb
temperature and state-wide energy consumption for the 2018 calendar year (ISO
New England, 2018), which we then aggregate into daily averages.

While the marginal distribution over daily average temperature is nearly
identical across states in the original dataset, the causal relationship
between temperature and energy consumption differs across states, likely due
to differences in population density, and commercial/industrial activity. To
introduce confounding, we systematically sample days (instances) from states
(objects) based on the state’s typical energy consumption, including hotter
days with higher probability for high consuming states. Specifically, we use
importance resampling with a target distribution over Farenheit temperatures
$T\sim\mathcal{N}(45+bias\cdot s_{o},15)$, where
$s_{CT}=3,s_{MA}=2,s_{ME}=1,s_{NH}=-1,s_{RI}=-2,s_{VT}=-3$. An example of this
sampling with $\text{bias}=9$ is shown in Figure 2 (a-c). Biased sampling in
this way introduces a statistical dependency across the dataset (consumption
is significantly higher in hotter days), that differs from the causal
relationship (low or high temperature causes a moderate increase in energy
consumption). This approach of sampling quasi-experimental data to simulate
confounding is an emerging standard in causal inference evaluation (Gentzel et
al., 2019) although existing benchmarks are not hierarchically structured.
Figure 2 (a-d) shows an example of this sampling process for the NEEC
benchmark.

Model | CT | MA | ME | NH | RI | VT
---|---|---|---|---|---|---
GP-SLC | 1.0 | 1.0 | 1.0 | 1.0 | 1.0 | 1.0
GP-NoConf | 13.2 | 13 | 31.5 | 41.6 | 47.4 | 14.9
GP-NoObj | 19.1 | 14 | 26.8 | 36.2 | 48 | 16.5
GP-PerObj | 1.6 | 1.3 | 5.2 | 9.7 | 6.5 | 0.7
MLM1 | 6.9 | 5 | 25.0 | 5 | 5.1 | 0.7
MLM2 | 6.4 | 4.9 | 39.4 | 6.3 | 9.9 | 3.8
BART | 4.1 | 2.1 | 13.3 | 3.6 | 3.3 | 2.4
Table 3: $\sqrt{\textrm{MSE}}$ for the New England energy consumption
benchmark, with $\text{bias}=9^{\circ}F$ and 25 samples per state. Lower is
better. Scores are normalized by GP-SLC’s score for the same state.

Sampling in this way does not provide instance-level counterfactual outcomes.
Instead, we estimate the sample-average ground truth counterfactual outcome by
fitting a Gaussian process regression model for each state, using treatments
and outcomes from the entire calendar year.

Figure 2(e) shows the models’ performances with varying degree of confounding
and sample sizes, and Figure 3 shows the estimated and actual effect of
temperate on electric energy consumption for two of the six states. Despite
the induced confounding, GP-SLC consistently produces accurate estimates of
causal effect. The baselines that ignore confounding (GP-NoConf and GP-NoObj)
perform poorly as the degree of confounding increases, incorrectly attributing
sample-wide association as indicative of causal effect. The linear multi-level
models (MLM 1 and MLM 2) are not biased by confounding, but produce poor
estimates due to their restrictive parametric assumptions. The remaining two
baselines (GP-PerObj and BART) produce more accurate estimates than the other
four baselines, but still overfit.

###  6.4 Limitations

Despite the fact that GP-SLC produces state-of-the-art counterfactual
predictions on most of our synthetic and semisynthetic benchmarks, it tends to
underestimate the uncertainty in these estimates. In other words, the
posterior density on the ground-truth counterfactual is sometimes low, despite
the fact that the mean estimate is close to the ground-truth relative to the
baselines. We suspect that this is partially attributable to inaccuracies
resulting from our approximate inference procedure (Algorithms 2 and 3).
Alternative approximate inference schemes, such as using our current approach
as a rejuvenation move in a sequential Monte Carlo (SMC) algorithm (Doucet et
al., 2001), may resolve these inaccuracies. This kind of SMC-based inference
procedure may also help GP-SLC scale to problems with more covariates and
objects than we explore in this paper.

Our empirical study focusses on data generating processes that satisfy GP-
SLC’s implicit semiparametric assumptions; (i) covariates for individuals
belonging to the same object are marginally Gaussian distributed, and (ii)
exogenous noise is additive and Gaussian. The effect of these modeling
assumptions on counterfactual prediction and estimates of effect strength
needs additional empirical characterization, ideally via large-scale synthetic
experiments (where ground truth is known and robustness to modeling bias can
be qualitatively studied).

##  7 Related Work

Leveraging hierarchical structure is well-established as a technique for
adjusting for latent confounding (Gelman, 2006; Gelman & Hill, 2006; Hong &
Raudenbush, 2006). Using Gaussian processes for causal inference is also well-
established (Alaa & van der Schaar, 2017, 2018; Silva & Gramacy, 2010; Schulam
& Saria, 2017; Zhang et al., 2010), as is the use of generative model
approaches to adjust for latent confounders given restrictions on structure
(Miao et al., 2018; Louizos et al., 2017; Tran & Blei, 2018; Wang & Blei,
2019). To the best of our knowledge, GP-SLC is the first semiparametric
generative modeling approach that leverages hierarchical structure to adjust
for latent confounders.

GP-SLC is one of many recent techniques (Shalit et al., 2017; Johansson et
al., 2016) for estimating individual-level treatment effects. Prior work
focusses on the propositional setting under strong ignorability, i.e. with no
latent confounders. We focus on the hierarchical setting in which latent
confounders are shared across multiple instances.

Recent work (Schulam & Saria, 2017) has used Gaussian process models for
causal inference in temporal settings, which assumes unconfoundedness and that
the outcome is smooth with respect to time and covariates. GP-SLC allows for
the existence of object-level latent confounders, and instead assumes that the
outcome is smooth with respect to treatment assignment, covariates, and latent
confounders. Longitudinal data analysis is closely related to the hierarchical
settings we consider in this work: measurements (instances) of individuals
(objects) are repeated over a period of time. Extending GP-SLC to the setting
where latent confounders are not shared across instances, but instead change
over time, is an exciting area of future work.

GP-SLC is most similar to (Alaa & van der Schaar, 2017), in that their
approach also uses GP models to estimate individual treatment effects.
However, GP-SLC: (i) handles hierarchical latent confounders by first
performing inference over object-level latent variables; (ii) accounts for the
covariance between noise-free factual and counterfactual outcomes (see
$\Sigma_{12}$ and $\Sigma_{21}$ in the Supplementary materials); and (iii)
uses a Monte Carlo algorithm for inference that yields quantified uncertainty
estimates. Their approach could be applied in hierarchical settings by
treating the object identifier $o$ as a categorical covariate and using a
delta kernel to construct the outcome kernel covariance matrix. This is
identical to the GP-PerObj baseline, except that GP-PerObj does not share
inferred kernel hyperparameters across objects.

##  8 Conclusions

This paper presents GP-SLC, a Gaussian process model for causal inference with
hierarchically structured latent confounders. In Section 6, we show that,
compared to widely used alternatives, GP-SLC produces more accurate estimates
of causal effect in realistic sparse observational settings where strong prior
knowledge about structure can inform causal estimates. The hierarchical
structure we exploit in this paper is one of many kinds of structural
background knowledge that could improve causal estimates, and developing
techniques to exploit such knowledge is an important area of future work.
Extending GP-SLC to handle large observational datasets (Cao, 2018; Quiñonero-
Candela & Rasmussen, 2005) or to leverage experimental evidence (Witty et al.,
2019) are also exciting areas of future work.

## Acknowledgments

Thanks to Marco Cusumano-Towner, Feras Saad, Alex Lew, Cameron Freer, Rachel
Paiste, Amanda Gentzel, Andy Zane, Jameson Quinn, and the anonymous reviewers
for their helpful feedback and suggestions. Sam Witty, Kenta Takatsu, and
David Jensen were supported by DARPA and the United States Air Force under the
XAI (Contract No. HR001120C0031) and CAML (Contract No. FA8750-17-C-0120)
programs, respectively. Vikash Mansinghka was supported by DARPA under the SD2
program (Contract No. FA8750-17-C-0239) and a philanthropic gift from the
Aphorism Foundation. Any opinions, findings and conclusions or recommendations
expressed in this material are those of the authors and do not necessarily
reflect the views of DARPA or the United States Air Force.

## References

  * Alaa & van der Schaar (2017) Alaa, A. and van der Schaar, M.  Bayesian inference of individualized treatment effects using multi-task Gaussian processes.  In _Advances in Neural Information Processing Systems_ , pp. 3424–3432, 2017.
  * Alaa & van der Schaar (2018) Alaa, A. and van der Schaar, M.  Bayesian nonparametric causal inference: Information rates and learning algorithms.  _IEEE Journal of Selected Topics in Signal Processing_ , 12(5):1031–1046, 2018.
  * Berkson (1946) Berkson, J.  Limitations of the application of fourfold table analysis to hospital data.  _Biometrics Bulletin_ , 2(3):47–53, 1946.
  * Boomsma et al. (2002) Boomsma, D., Busjahn, A., and Peltonen, L.  Classical twin studies and beyond.  _Nature Reviews Genetics_ , 3(11):872–882, 2002\.
  * Cao (2018) Cao, Y.  _Scaling Gaussian Processes_.  PhD thesis, University of Toronto (Canada), 2018.
  * Cusumano-Towner et al. (2019) Cusumano-Towner, M. F., Saad, F. A., Lew, A. K., and Mansinghka, V. K.  Gen: A general-purpose probabilistic programming system with programmable inference.  In _Proceedings of the 40th ACM SIGPLAN Conference on Programming Language Design and Implementation_ , pp. 221–236. ACM, 2019.
  * D’Amour (2019) D’Amour, A.  On multi-cause approaches to causal inference with unobserved counfounding: Two cautionary failure cases and a promising alternative.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pp. 3478–3486, 2019.
  * Doob (1949) Doob, J. L.  Application of the theory of martingales.  _Le calcul des probabilites et ses applications_ , pp. 23–27, 1949\.
  * Doucet et al. (2001) Doucet, A., De Freitas, N., and Gordon, N.  An introduction to sequential monte carlo methods.  In _Sequential Monte Carlo methods in practice_ , pp. 3–14. Springer, 2001.
  * Elwert & Winship (2014) Elwert, F. and Winship, C.  Endogenous selection bias: The problem of conditioning on a collider variable.  _Annual review of sociology_ , 40:31–53, 2014.
  * Gelman (2006) Gelman, A.  Multilevel (hierarchical) modeling: What it can and cannot do.  _Technometrics_ , 48(3):432–435, 2006.
  * Gelman & Hill (2006) Gelman, A. and Hill, J.  _Data Analysis Using Regression and Multilevel/Hierarchical Models_.  Cambridge University Press, 2006.
  * Gentzel et al. (2019) Gentzel, A., Garant, D., and Jensen, D.  The case for evaluating causal models using interventional measures and empirical data.  In _Advances in Neural Information Processing Systems_ , pp. 11717–11727, 2019.
  * Hastings (1970) Hastings, K.  _Monte Carlo Sampling Methods Using Markov Chains and Their Applications_.  Oxford University Press, 1970.
  * Hill (2011) Hill, J.  Bayesian nonparametric modeling for causal inference.  _Journal of Computational and Graphical Statistics_ , 20(1):217–240, 2011.
  * Hong & Raudenbush (2006) Hong, G. and Raudenbush, S.  Evaluating kindergarten retention policy.  _Journal of the American Statistical Association_ , 101(475):901–910, 2006.
  * Hong & Yu (2008) Hong, G. and Yu, B.  Effects of kindergarten retention on children’s social-emotional development: An application of propensity score method to multivariate, multilevel data.  _Developmental Psychology_ , 44(2):407, 2008.
  * Imbens & Rubin (2015) Imbens, G. and Rubin, D.  _Causal Inference in Statistics, Social, and Biomedical Sciences_.  Cambridge University Press, 2015.
  * ISO New England (2018) ISO New England.  Energy, load, and demand reports.  https://www.iso-ne.com/isoexpress/web/reports/load-and-demand/-/tree/zone-info, 2018\.
  * Jensen et al. (2019) Jensen, D., Burroni, J., and Rattigan, M.  Object conditioning for causal inference.  In _Proceedings of the Thirty-Fifth Conference on Uncertainty in Artificial Intelligence (UAI)_. AUAI Press, 2019.
  * Johansson et al. (2016) Johansson, F., Shalit, U., and Sontag, D.  Learning representations for counterfactual inference.  In _Proceedings of The 33rd International Conference on Machine Learning_ , volume 48 of _Proceedings of Machine Learning Research_ , pp. 3020–3029, New York, New York, USA, 20–22 Jun 2016. PMLR.
  * Kuroki & Pearl (2014) Kuroki, M. and Pearl, J.  Measurement bias and effect restoration in causal inference.  _Biometrika_ , 101(2):423–437, 2014.
  * Lawrence (2004) Lawrence, N.  Gaussian process latent variable models for visualisation of high dimensional data.  In _Advances in Neural Information Processing Systems_ , pp. 329–336, 2004.
  * Liang & Zeger (1986) Liang, K.-Y. and Zeger, S. L.  Longitudinal data analysis using generalized linear models.  _Biometrika_ , 73(1):13–22, 1986.
  * Loftus & Masson (1994) Loftus, G. and Masson, M.  Using confidence intervals in within-subject designs.  _Psychonomic Bulletin & Review_, 1(4):476–490, 1994.
  * Louizos et al. (2017) Louizos, C., Shalit, U., Mooij, J., Sontag, D., Zemel, R., and Welling, M.  Causal effect inference with deep latent-variable models.  In _Advances in Neural Information Processing Systems_ , pp. 6446–6456, 2017.
  * Miao et al. (2018) Miao, W., Geng, Z., and Tchetgen Tchetgen, E.  Identifying causal effects with proxy variables of an unmeasured confounder.  _Biometrika_ , 105(4):987–993, 08 2018.
  * Murray et al. (2010) Murray, I., Prescott Adams, R., and MacKay, D. J.  Elliptical slice sampling.  2010\.
  * Neal (2012) Neal, R.  _Bayesian Learning for Neural Networks_ , volume 118.  Springer Science & Business Media, 2012.
  * Pearl (2009) Pearl, J.  _Causality: Models, Reasoning and Inference_.  Cambridge University Press, New York, NY, USA, 2nd edition, 2009.
  * Pearl (2011) Pearl, J.  The algorithmization of counterfactuals.  _Annals of Mathematics and Artificial Intelligence_ , 61(1):29, 2011.
  * Quiñonero-Candela & Rasmussen (2005) Quiñonero-Candela, J. and Rasmussen, C. E.  A unifying view of sparse approximate gaussian process regression.  _Journal of Machine Learning Research_ , 6(Dec):1939–1959, 2005.
  * Ramey et al. (1992) Ramey, C. T., Bryant, D. M., Wasik, B. H., Sparling, J. J., Fendt, K. H., and La Vange, L. M.  Infant health and development program for low birth weight, premature infants: Program elements, family participation, and child intelligence.  _Pediatrics_ , 89(3):454–465, 1992.
  * Rasmussen (2003) Rasmussen, C.  Gaussian processes in machine learning.  In _Summer School on Machine Learning_ , pp. 63–71. Springer, 2003\.
  * Schulam & Saria (2017) Schulam, P. and Saria, S.  Reliable decision support using counterfactual models.  In _Advances in Neural Information Processing Systems_ , pp. 1697–1708, 2017.
  * Shadish et al. (2008) Shadish, W., Clark, M., and Steiner, P.  Can nonrandomized experiments yield accurate answers? a randomized experiment comparing random and nonrandom assignments.  _Journal of the American Statistical Association_ , 103(484):1334–1344, 2008.
  * Shalit et al. (2017) Shalit, U., Johansson, F., and Sontag, D.  Estimating individual treatment effect: Generalization bounds and algorithms.  In _Proceedings of the 34th International Conference on Machine Learning - Volume 70_ , ICML’17, pp. 3076–3085. JMLR.org, 2017.
  * Silva & Gramacy (2010) Silva, R. and Gramacy, R. B.  Gaussian process structural equation models with latent variables.  In _Proceedings of the Twenty-Sixth Conference on Uncertainty in Artificial Intelligence_ , pp. 537–545, 2010.
  * Titsias & Lawrence (2010) Titsias, M. and Lawrence, N.  Bayesian gaussian process latent variable model.  In _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics_ , pp. 844–851, 2010.
  * Tran & Blei (2018) Tran, D. and Blei, D. M.  Implicit causal models for genome-wide association studies.  In _International Conference on Learning Representations_ , 2018.
  * Van der Laan et al. (2007) Van der Laan, M. J., Polley, E. C., and Hubbard, A. E.  Super learner.  _Statistical applications in genetics and molecular biology_ , 6(1), 2007.
  * Wang & Blei (2019) Wang, Y. and Blei, D. M.  The blessings of multiple causes.  _Journal of the American Statistical Association_ , pp. 1–71, 2019\.
  * Witty et al. (2019) Witty, S., Lew, A., Jensen, D., and Mansinghka, V.  Bayesian causal inference via probabilistic program synthesis.  _arXiv preprint_ , arXiv:1910.14124, 2019.
  * Zhang et al. (2010) Zhang, K., Schölkopf, B., and Janzing, D.  Invariant Gaussian process latent variable models and application in causal discovery.  In _Proceedings of the 26th Conference on Uncertainty in Artificial Intelligence (UAI)_ , 2010.

##  9 Supplementary Materials

###  9.1 Kernel Functions

In this section, we present a detailed definition of each of the kernel
functions used in GP-SLC:

$\displaystyle
k^{\prime}_{x_{k}}([\textbf{u}_{o=Pa(i)}],[\textbf{u}_{o^{\prime}=Pa(i^{\prime})}])$
$\displaystyle=\sigma^{2}_{x}\exp{\left[-\sum_{j}^{N_{U}}\frac{(\textbf{u}_{o,j}-\textbf{u}_{o^{\prime},j})^{2}}{\lambda_{ux_{j,k}}}\right]}$
$\displaystyle
k^{\prime}_{t}([\textbf{u}_{o=Pa(i)},\mathrm{\textbf{x}}_{i}],[\textbf{u}_{o^{\prime}=Pa(i^{\prime})},\mathrm{\textbf{x}}_{i^{\prime}}])$
$\displaystyle=\sigma^{2}_{t}\exp{\left[-\sum_{j}^{N_{U}}\frac{(\textbf{u}_{o,j}-\textbf{u}_{o^{\prime},j})^{2}}{\lambda_{ut_{j}}}-\sum_{k}^{N_{X}}\frac{(\mathrm{\textbf{x}}_{i,k}-\mathrm{\textbf{x}}_{i^{\prime},k})^{2}}{\lambda_{xt_{k}}}\right]}$
$\displaystyle
k^{\prime}_{y}([\textbf{u}_{o=Pa(i)},\mathrm{\textbf{x}}_{i},\textbf{t}_{i}],[\textbf{u}_{o^{\prime}=Pa(i^{\prime})},\mathrm{\textbf{x}}_{i^{\prime}},\textbf{t}_{i^{\prime}}])$
$\displaystyle=\sigma^{2}_{y}\exp{\left[-\sum_{j}^{N_{U}}\frac{(\textbf{u}_{o,j}-\textbf{u}_{o^{\prime},j})^{2}}{\lambda_{uy_{j}}}-\sum_{k}^{N_{X}}\frac{(\mathrm{\textbf{x}}_{i,k}-\mathrm{\textbf{x}}_{i^{\prime},k})^{2}}{\lambda_{xy_{k}}}-\frac{(\textbf{t}_{i}-\textbf{t}_{i^{\prime}})^{2}}{\lambda_{ty}}\right]}.$

where $\lambda_{*}$ is a lengthscale hyperparameter and defined for each
dimension of corresponding variables. Here, each dimension of x is generated
independently given u, and $k^{\prime}_{x_{k}}$ refers to the kernel function
for the $k$th dimension of $x$. Intuitively, each kernel lengthscale
determines the relative strength of influence of each variable’s parents in
Equation 1. For example, if $\lambda_{ty}>>\lambda_{xy_{i=1...N_{X}}}$, the
covariance between instances (or counterfactuals) with similar treatments will
be greater than the covariance between instances with similar covariates.

###  9.2 Exact Inference:
$\boldsymbol{Y^{\prime}_{*}}-\boldsymbol{Y^{\prime}}$ Details

Here we provide additional details on how to compute GP-SLC’s conditional
distribution over individual treatment effects. Given the expression for
$\Bigg{(}\begin{bmatrix}[l]Y\\\ Y^{\prime}\\\
Y^{\prime}_{*}\end{bmatrix}\Bigg{|}\;T_{*},T,X,U,\Theta\Bigg{)}$ in Section
4.2, conditioning on $Y$ yields the following:

$\Big{(}\begin{bmatrix}[l]Y^{\prime}\\\
Y^{\prime}_{*}\end{bmatrix}\Big{|}T_{*},Y,T,X,U,\Theta\Big{)}\sim\mathcal{N}\Big{(}\begin{bmatrix}\mu_{1}\\\
\mu_{2}\end{bmatrix},\begin{bmatrix}\Sigma_{1,1}&\Sigma_{1,2}\\\
\Sigma_{2,1}&\Sigma_{2,2}\end{bmatrix}\Big{)}$

where,

$\displaystyle\mu_{1}$
$\displaystyle=K^{\prime}(\textbf{w},\textbf{w})K(\textbf{w},\textbf{w})\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}Y$
$\displaystyle\mu_{2}$
$\displaystyle=K^{\prime}(\textbf{w},\textbf{w}_{*})K(\textbf{w},\textbf{w})\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}Y$
$\displaystyle\Sigma_{1,1}$
$\displaystyle=K^{\prime}(\textbf{w},\textbf{w})-K^{\prime}(\textbf{w},\textbf{w})K(\textbf{w},\textbf{w})\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}K^{\prime}(\textbf{w},\textbf{w})$
$\displaystyle\Sigma_{1,2}$
$\displaystyle=K^{\prime}(W,W_{*})-K^{\prime}(W,W)K(W,W)\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}K^{\prime}(W,W_{*})$
$\displaystyle\Sigma_{2,1}$
$\displaystyle=K^{\prime}(W_{*},W)-K^{\prime}(W_{*},W)K(W,W)\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}K^{\prime}(W,W)$
$\displaystyle\Sigma_{2,2}$
$\displaystyle=K^{\prime}(W_{*},W_{*})-K^{\prime}(W_{*},W)K(W,W)\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}K(W,W_{*})$

As the difference of variables that are jointly Gaussian is Gaussian, we have
that
$(Y^{\prime}_{*}-Y^{\prime}|T_{*},X,T,Y,U,\Theta)\sim\mathcal{N}(\mu_{\textit{ITE}},\Sigma_{\textit{ITE}})$,
where $\mu_{\textit{ITE}}=\mu_{2}-\mu_{1}$ and
$\Sigma_{\textit{ITE}}=\Sigma_{1,1}-\Sigma_{1,2}-\Sigma_{2,1}+\Sigma_{2,2}$.

###  9.3 Asymptotic Posterior Consistency

Here we provide proofs for Proposition 5.1 and Theorems 5.2 and 5.3. The
analysis in this section follows the setup presented in (D’Amour, 2019), with
the inclusion of shared latent confounding amongst individual instances. We
omit covariates $X$ from this analysis and assume that $N_{U}=1$ for brevity
without loss of generality. Note that these theoretical results also hold for
the random intercepts multilevel model (Gelman, 2006).

Assuming linear kernels and additive Gaussian exogenous noise, we can
equivalently rewrite the GP-SLC model as follows. This equivalent structural
causal model is parameterized by latent variables
$\alpha,\beta,\tau\in\mathbb{R}$ and
$\sigma^{2}_{U},\sigma^{2}_{T},\sigma^{2}_{Y}\in\mathbb{R^{+}}$. For all $o\in
1,...,N_{O}$ and $i\in 1,...,N_{I}$, we have that:

$\begin{aligned}
\boldsymbol{\epsilon}_{u_{o}}&\sim\mathcal{N}(0,\sigma^{2}_{U})\\\
\boldsymbol{\epsilon}_{t_{i}}&\sim\mathcal{N}(0,\sigma^{2}_{T})\\\
\boldsymbol{\epsilon}_{y_{i}}&\sim\mathcal{N}(0,\sigma^{2}_{Y})\end{aligned}\quad\begin{aligned}
\textbf{u}_{o}&=\boldsymbol{\epsilon}_{u_{o}}\\\
\textbf{t}_{i}&=\alpha\textbf{u}_{o=Pa(i)}+\boldsymbol{\epsilon}_{t_{i}}\\\
\textbf{y}_{i}&=\beta\textbf{t}_{i}+\tau\textbf{u}_{o=Pa(i)}+\boldsymbol{\epsilon}_{y_{i}}.\end{aligned}$

In this setting, estimating individual treatment effect reduces to estimating
$\beta$, as
$\textbf{y}_{i,\textbf{t}_{*}}-\textbf{y}_{i}=\beta(\textbf{t}_{*}-\textbf{t}_{i})$.
We make the following observations.

Proposition 5.1 When $N_{O}=N_{I}$, $\textit{ITE}_{t_{*}}$ is not
asymptotically consistent $\forall t_{*}\in\mathbb{R}$.

For a detailed proof of Proposition 5.1, see Proposition 1 in (D’Amour, 2019).
In summary, they show that given any set of latent parameters
$\Theta=(\alpha,\beta,\tau,\sigma^{2}_{U},\sigma^{2}_{T},\sigma^{2}_{Y})$,
there exists an alternative set of parameters $\Theta^{\prime}$ such that
$P(T,Y|\Theta)=P(T,Y|\Theta^{\prime})$ and $\beta\neq\beta^{\prime}$. In other
words, the structural causal model forms a linear system of equations that is
rank-deficient. The set of parameters that satisfy this condition construct an
ignorance region.

Extending their results to the Bayesian setting, we have that for any two sets
of parameters $\Theta$ and $\Theta^{\prime}$ on the same ignorance region, the
posterior odds ratio reduces to the prior odds ratio,
$\frac{P(\Theta|T,Y)}{P(\Theta^{\prime}|T,Y)}=\frac{P(\Theta)P(T,Y|\Theta)}{P(\Theta^{\prime})P(T,Y|\Theta^{\prime})}=\frac{P(\Theta)}{P(\Theta^{\prime})}$.
By definition, $\Theta$ is not asymptotically consistent, as the posterior
$P(\Theta|T,Y)$ depends on the prior $P(\Theta)$. The problem of asymptotic
consistency can be mitigated when $N_{O}<N_{I}$.

Theorem 5.2 Assume there exists an object $o$ that is the parent of $n$
instances, $I^{\prime}=\\{i^{\prime}_{1},...,i^{\prime}_{n}\\}$. Then
$\textit{ITE}_{t_{*}}$ is asymptotically consistent as $n$ approaches
$\infty,\forall t_{*}\in\mathbb{R}$.

###### Proof.

For all $i^{\prime}\in I^{\prime}$, we have that
$\textbf{y}_{i^{\prime}}=\beta\textbf{t}_{i^{\prime}}+C+\boldsymbol{\epsilon}_{y_{i^{\prime}}}$
for some constant $C\in\mathbb{R}$. Therefore, the covariance between $T$ and
$Y$ in $I^{\prime}$ is uniquely given by $\beta$, i.e.
$cov(\textbf{t}_{i^{\prime}\in I^{\prime}},\textbf{y}_{i^{\prime}\in
I^{\prime}})=\beta$. Estimating the covariance of a bivariate normal has a
unique maximum likelihood solution. Therefore, by the Bernstein-von Mises
Theorem (Doob, 1949) we have that the posterior over $\beta$, and thus
$\textit{ITE}_{t_{*}}$, is asymptotically consistent as $n$ approach $\infty$.
∎

Theorem 5.3 Assume there exists $n$ objects
$\mathbb{O}=\\{o_{1},...,o_{n}\\}$, each of which are the unique parents of
$k\geq 2$ instances
$I^{\prime}_{o}=\\{i^{\prime}_{o,1},...,i^{\prime}_{o,k_{o}}\\}$. Then
$\textit{ITE}_{t_{*}}$ is asymptotically consistent as $n$ approaches
$\infty$.

###### Proof.

For all $o\in\mathbb{O}$, $j\in\\{1,...,k_{o}\\}$ let
$\textbf{t}^{\prime}_{i^{\prime}_{o,j}}=\textbf{t}_{i^{\prime}_{o,j}}-\bar{\textbf{t}}_{o}$
and
$\textbf{y}^{\prime}_{i^{\prime}_{o,j}}=\textbf{y}_{i^{\prime}_{o,j}}-\bar{\textbf{y}}_{o}$,
where $\bar{\textbf{t}}_{o}=\sum_{j}{\textbf{t}_{i^{\prime}_{o,j}}}/k_{o}$ and
$\bar{\textbf{y}}_{o}=\sum_{j}{\textbf{y}_{i^{\prime}_{o,j}}}/k_{o}$, i.e.,
the sample average over all instances that share a parent object. Therefore,
$\textbf{t}^{\prime}_{i^{\prime}_{o,j}}=\alpha\textbf{u}_{o}+\boldsymbol{\epsilon}_{t_{i^{\prime}_{o,j}}}-\sum_{j}(\alpha\textbf{u}_{o}+\boldsymbol{\epsilon}_{t_{i^{\prime}_{o,j}}})/k_{o}=\boldsymbol{\epsilon}_{t_{i^{\prime}_{o,j}}}-\sum_{j}{\boldsymbol{\epsilon}_{t_{i^{\prime}_{o,j}}}}/k_{o}$
and
$\textbf{y}^{\prime}_{i^{\prime}_{o,j}}=\beta(\alpha\textbf{u}_{o}+\boldsymbol{\epsilon}_{t_{i^{\prime}_{o,j}}})+\tau\textbf{u}_{o}+\boldsymbol{\epsilon}_{y_{i^{\prime}_{o,j}}}-\sum_{j}(\beta(\alpha\textbf{u}_{o}+\boldsymbol{\epsilon}_{t_{i^{\prime}_{o,j}}})+\tau\textbf{u}_{o}+\boldsymbol{\epsilon}_{y_{i^{\prime}_{o,j}}})/k_{o}=\beta\textbf{t}^{\prime}_{i^{\prime}_{o,j}}+\boldsymbol{\epsilon}_{y_{i^{\prime}_{o,j}}}-\sum_{i}{\boldsymbol{\epsilon}_{y_{i^{\prime}_{o,j}}}}/k_{o}$.
As $\boldsymbol{\epsilon}_{y_{i^{\prime}_{o,j}}}$ is independent of
$\textbf{t}^{\prime}_{i^{\prime}_{o,j}}$, we have that the covariance between
$\textbf{t}^{\prime}_{i^{\prime}_{o,j}}$ and
$\textbf{y}^{\prime}_{i^{\prime}_{o,j}}$ is equal to $\beta$. Therefore, the
problem of estimating $\beta$ reduces to estimating the covariance of a
bivariate normal distribution, $P(T^{\prime},Y^{\prime})$, which has a unique
maximum likelihood solution. As in the proof of Theorem 5.2, by the Bernstein-
von Mises Theorem (Doob, 1949) we have that the estimate of $\beta$, and thus
$\textit{ITE}_{t_{*}}$, is asymptotically consistent as $n$ approach $\infty$.
∎

###  9.4 Bayesian Linear Multilevel Model Baseline

One of the baselines we use in the experiments is Bayesian linear multilevel
models (Gelman, 2006). We implement two multilevel models, which introduce
varying degrees of shared parameters across objects. The first multilevel
model, also known as a random slope and intercepts model, (MLM 1) fits the
observations using the following structural equations.

$\displaystyle\sigma_{y}^{2}\sim\gamma\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}(\alpha_{\sigma_{y}},\beta_{\sigma_{y}})$
$\displaystyle\alpha\sim\mathcal{N}(\mu_{\alpha},\Sigma_{\alpha})$
$\displaystyle\beta_{o}\sim\mathcal{N}(\mu_{\beta},\sigma^{2}_{\beta})\text{
for $o=1\dots N_{O}$}$
$\displaystyle\eta_{o}\sim\mathcal{N}(\mu_{\eta},\sigma^{2}_{\eta})\text{ for
$o=1\dots N_{O}$}$
$\displaystyle\textbf{y}_{i}\sim\mathcal{N}(\beta_{o=Pa(i)}\textbf{t}_{i}+\alpha^{T}\mathrm{\textbf{x}}_{i}+\eta_{o=Pa(i)},\sigma_{y}^{2})$

This model allows varying intercepts $\eta$ and treatment effect $\beta$
across objects while assuming $\alpha$ is held constant across objects.

The second multilevel model, also known as the random intercepts model, (MLM
2) fits the observations using the following structural equations.

$\displaystyle\sigma_{y}^{2}\sim\gamma\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}(\alpha_{\sigma_{y}},\beta_{\sigma_{y}})$
$\displaystyle\alpha\sim\mathcal{N}(\mu_{\alpha},\Sigma_{\alpha})$
$\displaystyle\beta\sim\mathcal{N}(\mu_{\beta},\sigma^{2}_{\beta})$
$\displaystyle\eta_{o}\sim\mathcal{N}(\mu_{\eta},\sigma^{2}_{\eta})\text{ for
$o=1\dots N_{O}$}$
$\displaystyle\textbf{y}_{i}\sim\mathcal{N}(\beta\textbf{t}_{i}+\alpha^{T}\mathrm{\textbf{x}}_{i}+\eta_{o=Pa(i)},\sigma_{y}^{2})$

This model allows varying intercepts $\eta$ across objects while assuming
$\alpha$ and $\beta$ are held constant across objects. We implement both
models in Gen (Cusumano-Towner et al., 2019). For both models, we use
$\alpha_{\sigma_{y}}=4.0,\beta_{\sigma_{y}}=4.0,\mu_{(\cdot)}=0,\sigma^{2}_{\alpha}=3.0,\sigma^{2}_{\beta}=1.0,$
and $\sigma^{2}_{\eta}=10.0$ as priors.

###  9.5 Synthetic Experiments

We examine the finite-sample behavior of the GP-SLC model using two synthetic
datasets that match GP-SLC’s assumptions about the existence of object-level
latent confounders ($U$) that simultaneously influence instance-level observed
treatments (T), covariates ($X$), and outcomes ($Y$). The following structural
equations summarize the data generating process:

$\displaystyle W_{j}$ $\displaystyle\sim\mathcal{N}(0,1I_{3})\;\text{ for
}j=1,2,3$ $\displaystyle\textbf{u}_{o}$
$\displaystyle\sim\mathcal{N}(0,0.5I_{3})\text{ for $o=1\dots N_{O}$}$
$\displaystyle\mathrm{\textbf{x}}_{i}$
$\displaystyle=W\cdot\textbf{u}_{o=pa(i)}+\boldsymbol{\epsilon}_{x_{i}}\text{
where }\boldsymbol{\epsilon}_{x_{i}}\sim\mathcal{N}(0,0.5I_{3})\text{ for
}i=1...N_{I}$ $\displaystyle\textbf{t}_{i}$
$\displaystyle=g_{t}(\mathrm{\textbf{x}}_{i},\textbf{u}_{o=pa(i)})+\boldsymbol{\epsilon}_{t_{i}}\text{
where }\boldsymbol{\epsilon}_{t_{i}}\sim\mathcal{N}(0,0.5)\text{ for
}i=1...N_{I}$ $\displaystyle\textbf{y}_{i}$
$\displaystyle=g_{y}(\textbf{t}_{i},\mathrm{\textbf{x}}_{i},\textbf{u}_{o=pa(i)})+\boldsymbol{\epsilon}_{y_{i}}\text{
where }\boldsymbol{\epsilon}_{y_{i}}\sim\mathcal{N}(0,0.5)\text{ for
}i=1...N_{I}$

First, we draw u from a multivariate Gaussian distribution. Then, we generate
covariates x as linear combinations of u with additive exogenous noise. We
generate treatments t as a function ($g_{t}$) of x and u with additive noise.
Finally, we generate outcome y as a function ($g_{y}$) of x, t, and u with
additive noise. For multi-dimensional variables, x and u, we first apply the
nonlinear function to each dimension of x and u, then we aggregate them by
summing across dimensions.

The nonlinear treatment and outcome functions are shown in Table 4.

Dataset | $g_{t}(\mathrm{\textbf{x}},\textbf{u})$ | $g_{y}(\textbf{t},\mathrm{\textbf{x}},\textbf{u})$
---|---|---
Additive | $\sum_{j}\mathrm{\textbf{x}}_{*,j}\;\text{sin}(\mathrm{\textbf{x}}_{*,j})-\sum_{j}\textbf{u}_{*,j}\;\text{sin}(\textbf{u}_{*,j})$ | $\textbf{t}\text{sin}(2\textbf{t})+\sum_{j}\mathrm{\textbf{x}}_{*,j}\;\text{sin}(\mathrm{\textbf{x}}_{*,j})+3\sum_{j}\textbf{u}_{*,j}\;\text{sin}(\textbf{u}_{*,j})$
Multiplicative | $\frac{1}{10}(\sum_{j}\mathrm{\textbf{x}}_{*,j}\;\text{sin}(\mathrm{\textbf{x}}_{*,j}))(\sum_{j}\textbf{u}_{*,j}\;\text{sin}(\textbf{u}_{*,j}))$ | $\frac{1}{10}(\textbf{t}\text{sin}(2\textbf{t}))(\sum_{j}\mathrm{\textbf{x}}_{*,j}\;\text{sin}(\mathrm{\textbf{x}}_{*,j}))(\sum_{j}\textbf{u}_{*,j}\;\text{sin}(\textbf{u}_{*,j}))$
Table 4: The functional form of $T$ and $Y$ for 2 synthetic datasets with
continuous treatments and nonlinear outcome functions.