BAYESIMP: Uncertainty Quantification for Causal Data Fusion
Year: 2022

Abstract:
While causal models are becoming one of the mainstays of machine learning, the problem of uncertainty quantification in causal inference remains challenging. In this paper, we study the causal data fusion problem, where datasets pertaining to multiple causal graphs are combined to estimate the average treatment effect of a target variable. As data arises from multiple sources and can vary in quality and quantity, principled uncertainty quantification becomes essential. To that end, we introduce Bayesian Interventional Mean Processes, a framework which combines ideas from probabilistic integration and kernel mean embeddings to represent interventional distributions in the reproducing kernel Hilbert space, while taking into account the uncertainty within each causal graph. To demonstrate the utility of our uncertainty estimation, we apply our method to the Causal Bayesian Optimisation task and show improvements over state-of-the-art methods.

heading: 1 Introduction
text: Causal inference has seen a significant surge of research interest in areas such as healthcare [1], ecology [2], and optimisation [3]. However, data fusion, the problem of merging information from multiple data sources, has received limited attention in the context of causal modelling, yet presents significant potential benefits for practical situations [4, 5]. In this work, we consider a causal data fusion problem where two causal graphs are combined for the purposes of inference of a target variable (see Fig.1). In particular, our goal is to quantify the uncertainty under such a setup and determine the level of confidence in our treatment effect estimates.
Let us consider the motivating example in Fig.1, where a medical practitioner is investigating how prostate cancer volume is affected by a statin drug dosage. We consider the case where the doctor
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
∗Denotes equal contribution with alphabetical ordering
only has access to two separate medical studies describing the quantities of interest. On one hand we have observational data, from one medical study D1 [1], describing the causal relationship between statin level and prostate specific antigen (PSA), and on the other hand we have observational data, from a second study D2 [7], that looked into the link between PSA level and prostate cancer volume. The goal is to model the interventional effect between our target variable (cancer volume) and the treatment variable (statin). This problem setting is different from the standard observational scenario as it comes with the following challenges:
• Unmatched data: Our goal is to estimate E[cancer volume|do(statin)] but the observed cancer volume is not paired with statin dosage. Instead, they are related via a mediating variable PSA.
• Uncertainty quantification: The two studies may be of different data quantity/quality. Furthermore, a covariate shift in the mediating variable, i.e. a difference between its distributions in two datasets, may cause inaccurate extrapolation. Hence, we need to account for uncertainty in both datasets.
Formally, let X be the treatment (Statin), Y be the mediating variable (PSA) and T our target (cancer volume), and our aim is to estimate E[T |do(X)]. The problem of unmatched data in a similar context has been previously considered by [5] using a two-staged regression approach (X → Y and Y → T ). However, uncertainty quantification, despite being essential if our estimates of interventional effects will guide decision-making, has not been previously explored. In particular, it is crucial to quantify the uncertainty in both stages as this takes into account the lack of data in specific parts of the space. Given that we are using different datasets for each stage, there are also two sources of epistemic uncertainties (due to lack of data) as well as two sources of aleatoric uncertainties (due to inherent randomness in Y and T ) [8] . It is thus natural to consider regression models based on Gaussian Processes (GP) [9], as they are able to model both types of uncertainties. However, as GPs, or any other standard regression models, are designed to model conditional expectations only and will fail to capture the underlying distributions of interest (e.g. if there is multimodality in Y as discussed in [10]). This is undesirable since, as we will see, interventional effect estimation requires accurate estimates of distributions. While one could in principle resort to density estimation methods, this becomes challenging since we typically deal with a number of conditional/ interventional densities.
In this paper, we introduce the framework of Bayesian Interventional Mean Processes (BAYESIMP) to circumvent the challenges in the causal data fusion setting described above. BAYESIMP considers kernel mean embeddings [11] for representing distributions in a reproducing kernel Hilbert space (RKHS), in which the whole arsenal of kernel methods can be extended to probabilistic inference (e.g. kernel Bayes rule [12], hypothesis testing [13], distribution regression [14]). Specifically, BAYESIMP uses kernel mean embeddings to represent the interventional distributions and to analytically marginalise out Y , hence accounting for aleatoric uncertainties. Further, BAYESIMP uses GPs to estimate the required kernel mean embeddings from data in a Bayesian manner, which allows to quantify the epistemic uncertainties when representing the interventional distributions. To illustrate the quality of our uncertainty estimates, we apply BAYESIMP to Causal Bayesian Optimisation [15], an efficient heuristic to optimise objective functions of the form x∗ = arg minx∈X E[T |do(X) = x]. Our contributions are summarised below:
1. We propose a novel Bayesian Learning of Conditional Mean Embedding (BAYESCME) that allows us to estimate conditional mean embeddings in a Bayesian framework.
2. Using BAYESCME, we propose a novel Bayesian Interventional Mean Process (BAYESIMP) that allows us to model interventional effect across causal graphs without explicit density estimation, while obtaining uncertainty estimates for E[T |do(X) = x].
3. We apply BAYESIMP to Causal Bayesian Optimisation, a problem introduced in [15] and show significant improvements over existing state-of-the-art methods.
Note that [16] also considered a causal fusion problem but with a different objective. They focused on extrapolating experimental findings across treatment domains, i.e. inferring E[Y |do(X)] when only data from p(Y |do(S)) is observed, where S is some other treatment variable. In contrast, we focus on modelling combined causal graphs, with a strong emphasis on uncertainty quantification. While [17] considered mapping interventional distributions in the RKHS to model quantities such as E[T |do(X)], they only considered a frequentist approach, which does not account for epistemic uncertainties.
Notations. We denote X,Y, Z as random variables taking values in the non-empty sets X ,Y and Z respectively. Let kx : X × X → R be positive definite kernels on X with an associated RKHS Hkx . The corresponding canonical feature map kx(x′, ·) is denoted as φx(x′). Analogously for Y and Z. In the simplest setting, we observe i.i.d samples D1 = {xi, yi, zi}Ni=1 from joint distribution PXY Z which we con-
catenate into vectors x := [x1, ..., xN ]>. Similarly for y and z. For this work, X is referred as treatment variable, Y as mediating variable and Z as adjustment variables accounting for confounding effects. With an abuse of notation, features matrices are defined by stacking feature maps along the columns, i.e Φx := [φx(x1), ..., φx(xN )]. We denote the Gram matrix as Kxx := Φ>x Φx and the vector of evaluations kxx as [kx(x, x1), ..., kx(x, xN )]. We define Φy,Φz analogously for y and z.
Lastly, we denote T = f(Y ) + as our target variable, which is modelled as some noisy evaluation of a function f : Y → T on Y while being some random noise. For our problem setup we observe a second dataset of i.i.d realisations D2 = {ỹj , tj}Mj=1 from the joint PY T independent of D1. Again, we define ỹ := [ỹ1, ..., ỹM ]> and t := [t1, ..., tM ]> just like for D1. See Fig.2 for illustration.

heading: 2 Background
text: Representing interventional distributions in an RKHS1 has been explored in different contexts [18, 17, 19]. In particular, when the treatment is continuous, [17] introduced the Interventional Mean Embeddings (IMEs) to model densities in an RKHS by utilising smoothness across treatments. Given that IME is an important building block to our contribution, we give it a detailed review by first introducing the key concepts of do-calculus [20] and conditional mean embeddings [21].
2.1 Interventional distribution and do-calculus
In this work, we consider the structural causal model [20] (SCM) framework, where a causal directed acyclic graph (DAG) G is given and encodes knowledge of existing causal mechanisms amongst the variables in terms of conditional independencies. Given random variables X and Y , a central question in interventional inference [20] is to estimate the distribution p(Y |do(X) = x), where {do(X) = x} represents an intervention on X whose value is set to x. Note that this quantity is not directly observed given that we are usually only given observational data, i.e, data sampled from the conditional p(Y |X) but not from the interventional density p(Y |do(X)). However, Pearl [20] developed do-calculus which allows us to estimate interventional distributions from purely observational distributions under the identifiability assumption. Here we present the backdoor and front-door adjustments, which are the fundamental components of DAG based causal inference.
The backdoor adjustment is applicable when there are observed confounding variables Z between the cause X and the effect Y (see Fig. 3 (Top)). In order to correct for this confounding bias we can use the following equation, adjusting for Z as p(Y |do(X) = x) = ∫ Z p(Y |X = x, z)p(z)dz.
The front-door adjustment applies to cases when confounders are unobserved (see Fig. 3 (Bottom)). Given a set of front-door adjustment variables Z, we can again correct the estimate for the causal effect fromX to Y with p(Y |do(X) = x) = ∫ Z ∫ X p(Y |x
′, z)p(z|X = x)p(x′)dx′dz. We rewrite the above formulae in a more general form as we show below. For the remainder of the paper we will opt for this notation:
p(Y |do(X) = x) = EΩx [p(Y |Ωx)] = ∫ p(Y |Ωx)p(Ωx)dΩx (1)
For backdoor we have Ωx = {X = x, Z} and p(Ωx) = δxp(Z) where δx is the Dirac measure at X = x. For front-door, Ωx = {X ′, Z} and p(Ωx) = p(X ′)p(Z|X = x).
1We refer the reader to a detailed review of RKHS methods provided in the Appendix of [5]

heading: 2.2 Conditional Mean Embeddings
text: Kernel mean embeddings of distributions provide a powerful framework for representing probability distributions [11, 21] in an RKHS. In particular, we work with conditional mean embeddings (CMEs) in this paper. Given random variables X,Y with joint distribution PXY , the conditional mean embedding with respect to the conditional density p(Y |X = x), is defined as:
µY |X=x := EY |X=x[φy(Y )] = ∫ Y φy(y)p(y|X = x)dy (2)
CMEs allow us to represent the distribution p(Y |X = x) as an element µY |X=x in the RKHSHky without having to model the densities explicitly. Following [21], CMES can be associated with a Hilbert-Schmidt operator CY |X : Hkx → Hky , known as the conditional mean embedding operator, which satisfies µY |X=x = CY |Xφx(x) where CY |X := CY XC−1XX with CY X := EY,X [φy(Y ) ⊗ φx(X)] and CXX := EX,X [φx(X)⊗ φx(X)] being the covariance operators. As a result, the finite sample estimator of CY |X based on the dataset {x,y} can be written as:
ĈY |X = Φy(Kxx + λI)−1ΦTx (3)
where λ > 0 is a regularization parameter. Note that from Eq.3, [22] showed that the CME can be interpret as a vector-valued kernel ridge regressor (V-KRR) i.e. φx(x) is regressed to an element in Hky . This is crucial as CMEs allow us to turn the integration, in Eq.2, into a regression task and hence remove the need for explicit density estimation. This insight is important as it allows us to derive analytic forms for our algorithms. Furthermore, the regression formalism of CMEs motivated us to derive a Bayesian version of CME using vector-valued Gaussian Processes (V-GP), see Sec.3.

heading: 2.3 Interventional Mean Embeddings
text: Interventional Mean Embeddings (IME) [17] combine the above ideas to represent interventional distributions in RKHSs. We derive the front-door adjustment embedding here but the backdoor adjustment follows analogously. Denote µY |do(X)=x as the IME corresponding to the interventional distribution p(Y |do(X) = x), which can be written as:
µY |do(X)=x := ∫ Y φy(y)p(y|do(X) = x)dy = ∫ X ∫ Z (∫ Y φy(y)p(y|x′, z)dy ) ︸ ︷︷ ︸
CME µY |X=x,Z=z
p(z|x)p(x′)dzdx′
using the front-door formula with adjustment variable Z, and rearranging the integrals. By definition of CME ∫ φy(y)p(y|x′, z)dy = CY |X,Z(φx(x′)⊗ φz(z)) and linearity of integration, we have
= CY |X,Z (∫ X φx(x
′)p(x′)dx′︸ ︷︷ ︸ =µX
⊗ ∫ Z φz(z)p(z|x)dz︸ ︷︷ ︸
=µZ|X=x
) = CY |X,Z ( µX ⊗ µZ|X=x )
Using notations from Sec.2.1, embedding interventional distributions into an RKHS is as follows. Proposition 1. Given an identifiable do-density of the form p(Y |do(X) = x) = EΩx [p(Y |Ωx)], the general form of the empirical interventional mean embedding is given by,
µ̂Y |do(X)=x = ΦY (KΩx + λI) −1ΦΩx(x) > (4)
where KΩx = KXX KZZ and ΦΩx(x) is derived depending on p(Ωx). In particular, for backdoor adjustments, Φ(bd)Ωx (x) = Φ > XkX(x, ·) Φ>Z µ̂z and for front-door Φ (fd) Ωx (x) = Φ>X µ̂X Φ>Z µ̂Z|X=x.

heading: 3 Our Proposed Method
text: Two-staged Causal Learning. Given two independent datasets D1 = {(xi, zi, yi)}Ni=1 and D2 = {(ỹj , tj)}Mj=1, our goal is to model the average treatment effect in T when intervening on variable X , i.e model g(x) = E[T |do(X) = x]. Note that the target variable T and the treatment variable X are never jointly observed. Rather, they are linked via a mediating variable Y observed in both datasets.
In our problem setting, we make the following two assumptions: (A1) The treatment only affects the target through the mediating variable, i.e T ⊥⊥ do(X)|Y → P (T |do(X), Y ) = p(T |Y ), in other words, that in the true data generating model P(X, Y, Z, T), all causal paths from X to T are mediated through Y. and (A2) Function f given by f(y) = E[T |Y = y] belongs to an RKHSHky .2
We can thus express the average treatment effect as: g(x) = E[T |do(X) = x] = ∫ Y
E[T |do(X) = x, Y = y]︸ ︷︷ ︸ =E[T |Y=y], since T⊥⊥do(X)|Y p(y|do(X) = x)dy (5)
= ∫ Y f(y)p(y|do(X) = x)dy = 〈f, µY |do(X)=x〉Hky . (6)
The final expression decomposes the problem of estimating g into that of estimating the IME µY |do(X) (which can be done using D1) and that of estimating the integrand f : Y → T (which can be done using D2). Each of these two components can either be estimated using a GP or KRR approach (See Table 1). Furthermore, the reformulation as an RKHS inner product is crucial, as it circumvents the need for density estimation as well as the need for subsequent integration in Eq.6. Rather, the main parts of the task can now be viewed as two instances of regression (recall that mean embeddings can be viewed as vector-valued regression).
To model g and quantify its uncertainty, we propose 3 GP-based approaches. While the first 2 methods, Interventional Mean Process (IMP) and Bayesian Interventional Mean Embedding (BAYESIME) are novel derivations that allow us to quantify uncertainty from either one of the datasets, we treat them as intermediate yet necessary steps to derive our main algorithm, Bayesian Interventional Mean Process (BAYESIMP), which allows us to quantify uncertainty from both sources in a principled way. For a summary of the methods, see Fig.4 and Table 1. All derivations are included in the appendix.
Interventional Mean Process: Firstly, we train f as a GP using D2 and model µY |do(X)=x as VKRR using D1. By drawing parallels to Bayesian quadrature [24] and conditional mean process introduced in [25], the integral of interest g(x) = ∫ f(y)p(y|do(X) = x)dy will be a GP indexed by the treatment variable X . We can then use the empirical embedding µ̂Y |do(X) learnt in D1 to obtain an analytic mean and covariance of g.
Bayesian Interventional Mean Embedding: Next, to account for the uncertainty from D1, we model f as a KRR and µY |do(X)=x using a V-GP. We introduce our novel Bayesian Learning of Conditional Mean Embeddings (BAYESCME), which uses a nuclear dominant kernel [26] construction, similar to [27], to ensure that the inner product 〈f, µY |do(X)=x〉 is well-defined. As the embedding is a GP, the resulting inner product is also a GP and hence takes into account the uncertainty in D1.(See Prop. 4). Bayesian Interventional Mean Process: Lastly, in order to account for uncertainties coming from both D1 and D2, we combine ideas from the above IMP and BAYESIME. We place GPs on both f and µY |do(X) and use their inner product to model E[T |do(X)]. Interestingly, the resulting uncertainty can be interpreted as the sum of uncertainties coming from IMP and BAYESIME with an additional interaction term (See Prop.5).

heading: 3.1 Interventional Mean Process
text: Firstly, we consider the case where f is modelled using a GP and µY |do(X)=x using a V-KRR. This allows us to take into account the uncertainty from D2 by modelling the relationship between Y and
2We note that this two-staged setup resembles Instrumental Variable [23] (IV) regression. However, our general setup allows the IV and the treatment to be confounded, which is not the case in standard IV regression.
T using in a GP. Drawing parallels to Bayesian quadrature [24] where integrating f with respect to a marginal measure results into a Gaussian random variable, we integrate f with respect to a conditional measure, thus resulting in a GP indexed by the conditioning variable. Note that [25] studied this GP in a non-causal setting, for a very specific downscaling problem. In this work, we extend their approach to model uncertainty in the causal setting. The resulting mean and covariance are then estimated analytically, i.e without integrals, using the empirical IME µ̂Y |do(X) learnt from D1, see Prop.2. Proposition 2 (IMP). Given dataset D1 = {(xi, yi, zi)}Ni=1 and D2 = {(ỹj , tj)}Mj=1, if f is the posterior GP learnt from D2, then g = ∫ f(y)p(y|do(X))dy is a GP GP(m1, κ1) defined on the treatment variable X with the following mean and covariance estimated using µ̂Y |do(X) ,
m1(x) = 〈µ̂Y |do(x),mf 〉Hky = ΦΩx(x) >(KΩx + λI) −1Kyỹ(Kỹỹ + λfI) −1t (7)
κ1(x, x ′) = µ̂>Y |do(x)µ̂Y |do(x′) − µ̂ > Y |do(x)Φỹ(Kỹỹ + λI) −1Φ>ỹ µ̂Y |do(x′) (8)
= ΦΩx(x) >(KΩx + λI) −1K̃yy(KΩx + λI) −1ΦΩx(x ′) (9)
where µ̂Y |do(x) = µ̂Y |do(X)=x,Kỹy = Φ>ỹ Φy, mf and K̃yy are the posterior mean function and covariance of f evaluated at y respectively. λ > 0 is the regularisation of the CME. λf > 0 is the noise term for GP f . Ωx is the set of variables as specified in Prop.1.
Summary: The posterior covariance between x and x′ in IMP can be interpreted as the similarity between their corresponding empirical IMES µ̂Y |do(X)=x and µ̂Y |do(X)=x′ weighted by the posterior covariance K̃yy, where the latter corresponds to the uncertainty when modelling f as a GP in D2. However, since f only considers uncertainty in D2, we need to develop a method that allows us to quantify uncertainty when learning the IME from D1. In the next section, we introduce a Bayesian version of CME, which then lead to BAYESIME, a remedy to this problem.

heading: 3.2 Bayesian Interventional Mean Embedding
text: To account for the uncertainty in D1 when estimating µY |do(X), we consider a GP model for CME, and later extend to the interventional embedding IME. We note that Bayesian formulation of CMEs has also been considered in [28], but with a specific focus on discrete target spaces.
Bayesian learning of conditional mean embeddings with V-GP. As mentioned in Sec.2, CMEs have a clear "feature-to-feature" regression perspective, i.e E[φy(Y )|X = x] is the result of regressing φy(Y ) onto φx(X). Hence, we consider a vector-valued GP construction to estimate the CME.
Let µgp(x, y) be a GP that models µY |X=x(y). Given that f ∈ Hky , for 〈f, µgp(x, ·)〉Hky to be well defined, we need to ensure µgp(x, ·) is also restricted to Hky for any fixed x. Consequently, we cannot define a GP(0, kx ⊗ ky) prior on µgp as usual, as draws from such prior will almost surely fall outside Hkx ⊗Hky [26]. Instead we define a prior over µgp ∼ GP(0, kx ⊗ ry), where ry is a nuclear dominant kernel [26] over ky, which ensures that samples paths of µgp live in Hkx ⊗Hky almost surely. In particular, we follow a similar construction as in [27] and model ry as ry(yi, yj) = ∫ ky(yi, u)ky(u, yj)ν(du) where ν is some finite measure on Y . Hence we can now setup a vector-valued regression inHky as follows:
φy(yi) = µgp(xi, ·) + λ 1 2 i (10)
where i ∼ GP(0, r) are independent noise functions. By taking the inner product with φy(y′) on both sides, we then obtain ky(yi, y′) = µgp(xi, y′) + λ 1 2 i(y
′). Hence, we can treat k(yi, yj) as noisy evaluations of µgp(xi, yj) and obtain the following posterior mean and covariance for µgp. Proposition 3 (BAYESCME). The posterior GP of µgp given observations {x,y} has the following mean and covariance:
mµ((x, y)) = kxx(Kxx + λI) −1KyyR −1 yyryy (11)
κµ((x, y), (x ′, y′)) = kxx′ry,y′ − kxx(Kxx + λI)−1kxx′ryyR−1yyryy′ (12)
In addition, the following marginal likelihood can be used for hyperparameter optimisation,
−N 2
( log |Kxx + λI|+ log |R| ) − 1
2 Tr
( (Kxx + λI) −1KyyR −1 yyKyy ) (13)
Note that in practice we fix the lengthscale of ky and ry when optimising the above likelihood. This is to avoid trivial solutions for the vector-valued regression problem as discussed in [10]. The Bayesian version of the IME is derived analogously and we refer the reader to appendix due to limited space.
Finally, with V-GPS on embeddings defined, we can model g(x) as 〈f, µgp(x, ·)〉Hky , which due to the linearity of the inner product, is itself a GP. Here, we first considered the case where f is a KRR learnt from D2 and call the model BAYESIME. Proposition 4 (BAYESIME). Given dataset D1 = {(xi, yi, zi)}Ni=1 and D2 = {(ỹj , tj)}Mj=1, if f is a KRR learnt from D2 and µY |do(X) modelled as a V-GP using D1, then g = 〈f, µY |do(X)〉 ∼ GP(m2, κ2) where,
m2(x) = ΦΩx(x) >(KΩx + λI) −1KyyR −1 yyRyỹA (14)
κ2(x, x ′) = BΦΩx(x) >ΦΩx(x)− CΦΩx(x)>(KΩx + λI)−1ΦΩx(x′) (15)
where A = (Kỹỹ + λfI)−1t, B = A>RỹỹA and C = A>RỹyR−1yyRyỹA
Summary: Constants B and C in κ2 can be interpreted as different estimation of ||f ||Hky , i.e the RKHS norm of f . As a result, for problems that are “harder” to learn in D2, i.e. corresponding to larger magnitude of ||f ||Hky , will result into larger values of B and C. Therefore the covariance κ2 can be interpreted as uncertainty in D1 scaled by the difficulty of the problem to learn in D2.

heading: 3.3 Bayesian Interventional Mean Process
text: To incorporate both uncertainties in D1 and D2, we combine ideas from IMP and BAYESIME to estimate g = 〈f, µY |do(X)〉 by placing GPs on both f and µY |do(X). Again as before, a nuclear dominant kernel ry was used to ensure the GP f is supported onHky . For ease of computation, we consider a finite dimensional approximation of the GPs f and µY |do(X) and estimate g as the RKHS inner product between them. In the following we collate y and ỹ into a single set of points ŷ, which can be seen as landmark points for the finite approximation [29]. We justify this in the Appendix.
Proposition 5 (BAYESIMP). Let f and µY |do(X) be GPs learnt as above. Denote f̃ and µ̃Y |do(X) as the finite dimensional approximation of f and µY |do(X) respectively. Then g̃ = 〈f̃ , µ̃Y |do(X)〉 has the following mean and covariance:
m3(x) = ExKyŷK −1 ŷŷRŷỹ(Rỹỹ + λfI) −1t (16)
κ3(x, x ′) = ExΘ > 1 R̃ŷŷΘ1E > x′︸ ︷︷ ︸
Uncertainty fromD1
+ Θ (a) 2 Fxx′ −Θ (b) 2 Gxx′︸ ︷︷ ︸
Uncertainty fromD2
+ Θ (a) 3 Fxx′ −Θ (b) 3 Gxx′︸ ︷︷ ︸
Uncertainty from Interaction
(17)
where Ex = ΦΩx(x) >(KΩx + λI) −1, Fxx′ = ΦΩx(x) >ΦΩx(x ′), Gxx′ = ΦΩx(x) >(KΩx + λI)−1ΦΩx(x ′), and Θ1 = K−1ŷŷRŷyR −1 yyKyy, Θ (a) 2 = Θ > 4 RŷŷΘ4,Θ (b) 2 = Θ > 4 RŷyR −1 yyRyŷΘ4 and Θ(a)3 = tr(K −1 ŷŷRŷŷK −1 ŷŷ R̄ŷŷ),Θ (b) 3 = tr(RŷyR −1 yyRyŷK −1 ŷŷ R̄ŷŷK −1 ŷŷ ) and Θ4 = K−1ŷŷRŷỹ(Kỹỹ + λf ) −1t. R̄ŷŷ is the posterior covariance of f evaluated at ŷ
Summary: While the first two terms in κ3 resemble the uncertainty estimates from IMP and BAYESIME, the last term acts as an extra interaction between the two uncertainties from D1 and D2. We note that unlike IMP and BAYESIME, g̃ from Prop.5 is not a GP as inner products between Gaussian vectors are not Gaussian. Nonetheless, the mean and covariance can be estimated.

heading: 4 Experiments
text: In this section, we first present an ablation studies on how our methods would perform under settings where we have missing data parts at different regions of the two datasets. We then demonstrate BAYESIMP’s proficiency in the Causal Bayesian Optimisation setting.
In particular, we compare our methods against the sampling approach considered in [15]. [15] start by modelling f : Y → T as GP and estimate the density p(Y |do(X)) using a GP along with do-calculus. Then given a treatment x, we obtain L samples of yl and R samples of fr from their posterior GPs.
The empirical mean and standard deviation of the samples {fr(yl)}L,Rl=1,r=1 can now be taken to estimate E[T |do(X) = x] as well as the correspondingly uncertainty. We emphasize that this point estimation requires repeated sampling and is thus inefficient compared to our approaches, where we explicitly model the uncertainty as covariance function.
Ablation study. In order to get a better intuition into our methods, we will start off with a preliminary example, where we investigate the uncertainty estimates in a toy case. We assume two simple causal graphs X −→ Y for D1 and Y −→ T for D2 and the goal is to estimate E[T |do(X) = x] (generating process given in the appendix). We compare our methods from Sec.3 with the sampling-based uncertainty estimation approach described above. In Fig.5 we plot the mean and the 95% credible interval of the resulting GP models for E[T |do(X) = x]. On the x-axis we also plotted a histogram of the treatment variable x to illustrate its density.
From Fig.5(a), we see that the uncertainty for sampling is rather uniform across the ranges of x despite the fact we have more data around x = 0. This is contrary to our methods, which show a reduction of uncertainty at high x density regions. In particular, x = −5 corresponds to an extrapolation of data, where x gets mapped to a region of y where there is no data in D2. This fact is nicely captured by the spike of credible interval in Fig.5(c) since IMP utilises uncertainty from D2 directly. Nonetheless, IMP failed to capture the uncertainty stemming from D1, as seen from the fact that the credible interval did not increase as we have less data in the region |x| > 5. In contrast, BAYESIME (Fig.5(d)) gives higher uncertainty around low x density regions but failed to capture the extrapolation phenomenon. Finally, BAYESIMP Fig.5(e) seems to inherit the desirable characteristics from both IMP and BAYESIME, due to taking into account uncertainties from both D1, D2. Hence, in the our experiments, we focus on BAYESIMP and refer the reader to the appendix for the remaining methods.
BayesIMP for Bayesian Optimisation (BO). We now demonstrate, on both synthetic and real-world data, the usefulness of the uncertainty estimates obtained using our methods in BO tasks. Our goal is to utilise the uncertainty estimates to direct the search for the optimal value of E[T |do(X) = x] by querying as few values of the treatment variable X as possible, i.e. we want to optimize for x∗ = arg minx∈X E[T |do(X) = x]. For the first synthetic experiment (see Fig.6 (Top)), we will use the following two datasets: D1 = {xi, ui, zi, yi}Ni=1 and D2 = {ỹj , tj}Mj=1. Note that BAYESIMP from Prop.5 is not a GP as inner products between Gaussian vectors are not Gaussian. Nonetheless, with the mean and covariance estimated, we will use moment matching to construct a GP out of BAYESIMP for posterior inference. At the start, we are given D1 and D2, where these observations are used to construct a GP prior for the interventional effect of X on T , i.e E[T |do(X) = x], to provide a “warm” start for the BO.
Again we compare BAYESIMP with the sampling-based estimation of E[T |do(X)] and its uncertainty, which is exactly the formulation used in the Causal Bayesian Optimisation algorithm (CBO) [15]. In order to demonstrate how BAYESIMP performs in the multimodal setting, we will be considering the case where we have the following distribution on Y i.e. p(y|u, z) = πp1(y|u, z) + (1− π)p2(y|u, z) where Y is a mixture and π ∈ [0, 1]. These scenarios might arise when there is an unobserved binary variable which induces a switching between two regimes on how Y depends on (U,Z). In this case, the GP model of [15] would only capture the conditional expecta-
tion of Y with an inflated variance, leading to slower convergence and higher variance in the estimates of the prior as we will show in our experiments. Throughout our experiments, similarly to [15], we will be using the expected improvement (EI) acquisition function to select the next point to query.
Synthetic data experiments. We compare BAYESIMP to CBO as well as to a simple GP with no learnt prior as baseline. We will be using N = 100 datapoints for D1 and M = 50 datapoints D2. We ran each method 10 times and plot the resulting standard deviation for each iteration in the figures below. The data generation and details were added in the Appendix. We see from the Fig.7 that
BAYESIMP is able to find the maxima much faster and with smaller standard deviations, than the current state-of-the-art method, CBO, using both front-door and backdoor adjustments (Fig.7(Right, Middle)). Given that our method uses more flexible representations of conditional distributions, we are able to circumvent the multimodality problem in Y . In addition, we also consider the unimodal version, i.e. π = 0 (see right Fig.7). We see that the performance of CBO improves in the unimodal setting, however BAYESIMP still converges faster than CBO even in this scenario.
Next, we consider a harder causal graph (see Fig.6 (Bottom)), previously considered in [15]. We again introduce multimodality in the Y variable in order to explore the case of more challenging conditional densities. We see from Fig.8 (Left, Middle), that BAYESIMP again converges much faster to the true optima than CBO [15] and the standard GP prior baseline. We note that the fast convergence of BAYESIMP throughout our experiments is not due to simplicity of the underlying BO problems. Indeed, the BO with a standard GP prior requires significantly more iterations. It is rather the availability of the observational data, allowing us to construct a more appropriate prior, which leads to a “warm” start of the BO procedure.
Healthcare experiments. We conclude with a healthcare dataset corresponding to our motivating medical example in Fig.1. The causal mechanism graph, also considered in the CBO paper [15], studies the effect of certain drugs (Aspirin/Statin) on Prostate-Specific Antigen (PSA) levels [6]. In our case, we modify statin to be continuous, in order to optimize for the correct drug dosage. However, in contrast to [15], we consider a second experimental dataset, arising from a different medical study, which looks into the connection between PSA levels and cancer volume amount in patients [7]. Similar to the original CBO paper [15], given that interventional data is hard to obtain, we construct data generators based on the true data collected in [7]. This is done by firstly fitting a GP on the data and then sampling from the posterior (see Appendix for more details). Hence this is the
perfect testbed for our model where we are interested in E[Cancer Volume|do(Statin)]. We see from Fig.8 (Right) that BAYESIMP again converges to the true optima faster than CBO hence allowing us to find the fastest ways of optimizng cancer volume by requesting much less interventional data. This could be critical as interventional data in real-life situations can be very expensive to obtain.

heading: 5 Discussion and Conclusion
text: In this paper we propose BAYESIMP for quantifying uncertainty in the setting of causal data fusion. In particular, our proposed method BAYESIMP allows us to represent interventional densities in the RKHS without explicit density estimation, while still accounting for epistemic and aleatoric uncertainties. We demonstrated the quality of the uncertainty estimates in a variety of Bayesian optimization experiments, in both synthetic and real-world healthcare datasets, and achieve significant improvement over current SOTA in terms of convergence speed. However, we emphasize that BAYESIMP is not designed to replace CBO but rather an alternative model for interventional effects.
In the future, we would like to improve BAYESIMP over several limitations. As in [15], we assumed full knowledge of the underlying causal graph, which might be limiting in practice. Furthermore, as the current formulation of BAYESIMP only allows combination of two causal graphs, we hope to generalise the algorithm into arbitrary number of graphs in the future. Causal graphs with recurrent structure will be an interesting direction to explore. Lastly, we would also like to include a cost function as in [15] to constrain the search space to the most sensible solutions.

heading: 6 Acknowledgements
text: The authors would like to thank Bobby He, Robert Hu, Kaspar Martens, Jake Fawkes and Joost van Amersfoort for helpful comments. SLC and JFT are supported by the EPSRC and MRC through the OxWaSP CDT programme EP/L016710/1. YWT and DS are supported in part by Tencent AI Lab and DS is supported in part by the Alan Turing Institute (EP/N510129/1). YWT’s research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) ERC grant agreement no. 617071.

Review summaries:
The paper proposes a Bayesian approach for estimating the treatment effect by combining two data sources. The main advantage of the Bayesian approach over the frequentist counterpart [17] is its ability to quantify uncertainty. The main machinery was GP which replaces KRR used in existing works. 

The paper studies the problem of uncertainty quantification for causal data fusion.  They introduce Bayesian Interventional Mean Processes, a framework that represents the interventional distributions using kernel mean embeddings and probabilistic integration. The framework is able to measure the uncertainty within the causal graph in each dataset. In experiments, the author demonstrates the performance improvement of the proposed method on the causal Bayesian optimization task.

The paper proposes a framework for embedding interventional distributions into RKHS that accounts for uncertainty in causal datasets especially when the confounders, treatments, mediators, and outcomes are not all observed in a single observational dataset. The main problem that is addressed here is quantifying the uncertainty in the estimation of causal effects when the treatment and the outcome are measured in two separate sources (potentially under different data-generating conditions) and the aim is to estimate the effect from a causal data fusion perspective. The paper introduces Bayesian Interventional Mean Processes for representing interventional distributions in the reproducing kernel Hilbert space. The method is applied to the task of Causal Bayesian Optimization where it is shown that the method converges faster as compared to other baseline methods. The main contribution of the paper is obtaining the uncertainty estimates for $E\[T\mid do(X) = x\]$.