Counterfactual Learning with Multioutput Deep Kernels
Authors: Alberto Caron
Ioanna Manolopoulou
Gianluca Baio
Venue: Accepted by TMLR ()

Abstract: In this paper, we address the challenge of performing counterfactual inference with observational data via Bayesian nonparametric regression adjustment, with a focus on high-dimensional settings featuring multiple actions and multiple correlated outcomes. We present a general class of counterfactual multi-task deep kernels models that estimate causal effects and learn policies proficiently thanks to their sample efficiency gains, while scaling well with high dimensions. In the first part of the work, we rely on Structural Causal Models (SCM) to formally introduce the setup and the problem of identifying counterfactual quantities under observed confounding. We then discuss the benefits of tackling the task of causal effects estimation via stacked coregionalized Gaussian Processes and Deep Kernels. Finally, we demonstrate the use of the proposed methods on simulated experiments that span individual causal effects estimation, off-policy evaluation and optimization.

Published in Transactions on Machine Learning Research (11/2022)
Counterfactual Learning with Multioutput Deep Kernels
Alberto Caron
alberto. cecertt.ccron... 19 19@ucll.cc. uk
Deportment of Statistical Science, University College London The Alan Turing Institute, London, UK
Gianluca Baio
g. babo wwwl. uk
Department of Statistioal Science, University College London
Ioanna Manolopoulou
i. Lmanolopouiouuuucc.ccco
Deportment of Statistical Science, University College London
Reviewed OpenReview: https: // opeenrve/pperevveee. net/ pet/forum? id= &GREAJdULX
Abstract
In this paper, we address the challenge of performing counterfactual inference with ob- servational data via Bayesian nonparametric regression adjustment, with a focus on high- dimensional settings featuring multiple actions and multiple correlated outcomes. We present a general class of counterfactual multi-task deep kernels models that estimate causal effects and learn policies proficiently thanks to their sample efficiency gains, while scaling well with high dimensions. In the first part of the work, we rely Structural Causal Models (SCM) to formally introduce the setup and the problem of identifying counterfactual quantities under observed confounding. We then discuss the benefits of tackling the task of causal effects estimation via stacked coregionalized Gaussian Processes and Deep Kernels. Finally, we demonstrate the use of the proposed methods on simulated experiments that span individual causal effects estimation, off-policy evaluation and optimization.
1 Introduction
In the recent years, there has been surge of attention towards the use of machine learning methods for causal inference under observational data. The desire for highly personalized crcision-making is indeed pervasive in many disciplines such as precision medicine, web advertising and education (Hodson 2016). However, in these fields, exploration of policies in the real-world is usually very costly and potentially harmful. For example, a randomized clinical trial targeted assessing the effects of a new surgical technique entails high expenses in terms of patient recruitment and design, in addition to concerns about safety and potential side effeets both treated and non-treated individuals. Thus, in the attempt to answer counterfactual questions about policy interventions, such "what would have happened if individual i undertook treatment. A instead of treatment B?', one cannot typically rely on randomized experimental data. On the other hand, observational data are abundant and more readily accessible at relatively lower costs, although they suffer from sample selection bias issues arising because of confounding factors. Counterfactual quantities can still be identified and estimated settings with observed (and in some cases unobserved) confounders, provided that some assumptions hold.
This work is aimed especially at tackling the problem of carrying out counterfactual inference using ob- servational data, with a focus on high-dimensional settings characterized by large number of covariates, multiple discrete actions (or manipulative variables) and multiple correlated outcomes of interest. To this end, we present a general class of counterfactual multi-task Deep Kernel Learning models (CounterDKL) that efficiently adapt to multiple actions and outcomes scttings by exploiting existing correlations, while inheriting the appealing scaling properties of DKL (Wilson et al. 2016) under large samples and large number of predictors, and preserving the Bayesian uncertainty quantification properties of Gaussian Processes (GPs).

---

Published in Transactions on Machine Learning Research (11/2022)
CounterDKL essentially consists of two joint, end-to-end, components: i) deep learning architecture that learns a lower-dinensional representation of high-dimensional input space: ii) multitask GP (kernel-based) component (Teh et al. 2005 Bonilla et al. 2008 Alvarez et al., 2012; Bohn et al, 2019) placed on the lower-dimensional representations that can learn the posterior joint distribution of the outcomes conditional on inpuls and actions, by avoiding parameter proliferation and stability issues that arise from placing a multitask GP kernel directly on the high-dimensional input space. Specifically, our contributions can be listed as follows:
We extend the class of causal multitask GPs (Alaa & van der Schaar 2017 2018) to multiple actions and multiple outcomes designs, by specifying stacked coregionalization model, highlighting its advantages in tackling "poor overlap' regions and disadvantages in terms of scalability.
We introduce the class of CounterDKL, for causal inference under multiple actions/ /outooossss discussing their benefits over causal multitask GPs in high-dimensional settings (in terms of covariates, actions and outcomes).
We demonstrate the use of CounterDKL, with simulated experiments on causal effects estimation, off- policy evaluation (OPE) and learning off-policy (OPL) problems (Dudík et al. 2011 Dudík et al. 2014
Farajtabar et al. 2018 Kallus 2021 by providing also an Python implement apppemmnntaiion ation of the models, based
on GPyTorch
Studies with multiple outcomes are quile common in applied research, as policy decisions are rarely based on a single outcome, but rather on a profile of different outcomes that might exhibit positive or negative correlation. As an example, the medical domain, prescription of a specific treatment, such as anti-coagulants to lower cholesterol level, depends on the risk of myocardial infarction (primary outcome), as well as the risk of bleeding (unwanted, correlated, side effect).
1.1 Related Work
Many significant contributions in the literature on non-linear regression techniques for countrfacttul//casst learning have been made in the recent years. These include works heterogeneous (or individual) treatment
effects estimation (Shalit et al. 2017
Künzel et al.
Yao et al. and on the related field of off-policy evaluation and learning
2017
Alaa & van der Schaar 2017 2018
2018 Wager& & Athey 2018 Nie el al.
2020
2018 Athey & Wagerr 2021), ,, that focus more on the prescriptive, rather than predictive, goal. Among the several contributions, we particularly draw attention to the early work of Hill (2011) that first proposed Bayesian nonparametric methods (specifically Bayesian Additive Regression Trees) for regression adjustment in estimating causal effects, followed by the extension of Hahn et al. (2020); Caron et al. (2022b), and the seminal work of Alaa & van der Schaar (2017 2018)
Dudik et al. 20144] Farajtabar et al.] 2018 Kallus
who first proposed multitask Gaussian Processes for causal effects estimation tasks, but limited to contexts with binary actions and a single (continuous) outcome. All these works highlight. the advantages of Bayesian nonparametric models in terms of flexible non-lincar function approximation and uncertainty quantification.
More specifically, the advantage of multitask GPs and DKL for causal learning lies in their sample efficiency gains. As observational data often feature imbalanced action arms with relatively few instances, sample splitting can result in under over-fitting of the surfaces of interest. Multitask GPs and DKL both allow for information ion sharing when learning actions and outcomes correlated tasks, while only DKL guarantees better scalability and also makes the choice of a GP kernel less challenging, as the deep learning structure can itself learn arbitrarily complex non-linear functions (Wilson et al. 2016
2 Problem Framework
In order to introduce the problems of causal effects identification and estimation, we borrow the notation from do-caleulus and Structural Causal Models (SCM) (Pearl. 2009), and start by defining the latter: Definition 2.1 (SCM). A Structural Causal Model is 4-tuple (E, V. F, F(E()) consisting of (subscript j indicates a random element in the set):
1 Full code at: http:::/ggttaab.ccm//aaaaaaaraea
2

---

Published in Transactions on Machine Learning Research (11/2022)
1. E: denoting a set of exogenous variables, defined variables determined outside of the model.
2. V: denoting a set of endogenous variables, defined as variables determined inside the model.
3. Fx a set of functions/ f, E F mapping each element E E E and every parent variables of V, E V which we denote by pa(V;) - to the endogenous variables V, E V, f3 Ej U pa(VV) - V.
4. p(e,): a probability distribution over E EU.
A SCM implies causal Directed Acyclic Graph (DAG), where nodes in the graph depict variables, i.c. Ej and v,, while edges denote the functional causal relationships f3 E F. In the following paragraphs, we examine the problem of identifying and estimating these causal relationships fy E F, under observational data, relying on SCMs and causal DAGs.
Suppose we have access to an observational dataset Dj = {X,, A,, Y, ~ p(-), with E NY, where X; E X is set of covariates, A, E A a set of discrete actions (or manipulative variables), and Y E R. a set of M different outcomes. Here we consider continuous type of outcomes for simplicity, but we extend the methods also to discrete outcomes (Milios et al. 2018); as in the experimental Sec- tion 5.2 The overarching goal is to identify and estimate the (average) effects of intervening on the manipulative variable A,, by setting it equal to some value a, on the outcomes Y, We denote the main quantity of interest, the joint interventional probability distribution of all the outcomes conditional on A, = a by using the do-operator as p(Y (do(A a)). We assume that the SCM is fully described by the following set of oquations:
X
Y
Causal DAG depicting the setting of interest, where arrows represent causal relation- ships. Set of observable covariales X E X satis- fies the backdoor criterion, in that conditioning on them allows for identification of the causal effects of action E A on the outcomes Y V.
Figure 1:
X -x=xx((,xx,
A = FA (p(( = =(,( Y fr (pa(Y)), (r,
(1)
where: (-) a vector-valued function if multiple outcomes are considerod; pa(A;) denotes parent variables (or causes) of Ar; Eij are error terms with a distribution p(εx,). The equivalent causal DAG is represented in Figure []] although we note that the methods presented in this work extend to other types of causal DAGs. We make two standard assumption for identification of the causal effect A + Y in this scenario. The first is unconfoundedness, stating that there are no unobserved confounders, equivalently that the set of observed covariates X; E X' is causally sufficient, in the sense that conditioning on X, allows to identify the causal association between A, and Y, Using Pearl's terminology, we equivalently say that X satisfics the backdoor criterion (see A.1 in the appendix for a formal definition), in that it "blocks all backdoor paths' from A to all the outcomes Y (Pearl 2009). The second assumption is overlap (also common support assumption). Overlap requires that 0 p(A; a|X; x)) 1, Va E - i.e. that the observed actions allocation given X; = is never deterministic, so we could theoretically observe data points for which X; = in each of the discrete arms of A (A.2 in the appendix). This ensures that we have comparable units in terms of X, in each action arm, so we can approximate (X, A, = a, Ei,)) well enough. Violation of overlap for portions of X' undermines generalization and extrapolation model's prediction those regions; thus, one must be careful as to which subpopulation, defined by a common support Xover C X, to target estimate causal effects. Under these two assumptions, the multivariate interventional distribution p( do(A o), x)) can be recovered via backdoor adjustment
known as positivity
p(Y do(A = a), X x) p(Y A = X x) (Pearl 2009) (see Appendix A). Hence, provided that the covariates X,, direct common causes of A, and Y, satisfy the backdoor criterion, we can estimate unbiased causal effects with the observed quantities in D, (X,, A,,
Extensions settings with contimuous action space A C R fentailing dose-response curve estimation) are non-trivial Tmail van Dyk 2004), and the multitask learning paradigm not specifically suitable for them.
3

---

Published in Transactions on Machine Learning Research (11/2022)
3
Counterfactual Learning with Multitask GPs
For ease of exposition, let us consider the simple case depicted by the causal DAG in Figure D with single
continuous outcome Y; E R, with i E [1,....N] We tackle the problem of estimating p doo = a)) via
non- non-linear
[Johansson et al. 2016 Shalit et al., 2017 Künzel et al. 2017 Nie & Wager
2020
Caron ct al. 2022a). In particular, we assumc, in line with most of the previous works, additive noisc such that the oulcome functional can be written as:
structure 3
V = (X; A.) + +,Y,, [e(eyy) = 00
There are different ways in which one can derive an estimator for do = a)) and its moments,
e.g. [(Y]do((4 =)), from 2 (Künzel et al. 2017 Caron et al. 2022a). Alaa & van der Schaar (2017
2018) first proposed the use of multitask learning via Gaussian Process regression, in the specific context.
of conditional average treatment effects estimation, which is defined, assuming binary A, E {0, 1), as the 0), æz... The idea behind causal multitask GPs is
quantity 7 r(x,)= do = =1), =)] [[
view the D AA interventional quantities Ya: where D the number of discrete action arms, the outpul. from vector-valued function fr() X + RD
(plus noise), modelled with GP prior:
fr()) ~ 100~GP(( SPPPPpPmp
with mean m(x;)= E R and covariance/kenn /ernel function K( 23) E RD x RD given two P-dimensional input points 2; 2j E X for units i and j. Given the likelihood function as a multivariate Gaussian p(y:lfr (i_ff,aa Fi, Σ) A (ff (, E), where E E TR. x RRD is the error covariance diagonal matrix with {o%) on the diagonal and E RD an output point, the posterior predictive distribution for a train set covariate realization I E X, train set outcome realization y: E R and a test set covariate realization xÿ E X is obtained as, assuming vero prior mean m(.) = 0 for simplicity:
p(fr (m)) (#: v.), fr, 0) = N K (2),
= x; (= K(tj, K(xj. x. (
(4)
where H = K(x;; xx) +
and where & denotes the model parameters and fy (xx)) the function evaluated at a test point æ). Under zero prior mean m(-) = 0, the multitask GP in 3 is fully parametrized by its kernel function K(;-)) The structure of the kernel function in multitask GP is what induces task-relatedness when fitting the multi-valued surface fy (-).
3.1 The multitask kernel
The simplest specification for the multitask kernel matrix is given by the separable kernels structure, which assumes single entries in ((.,)) to be of the form ka,a' (a,, (z) k(x;,, k.) = k(æ;, 2x,) ba,a', with action a E A = {l,..., D}. Here, k(w, æ;) represents a base kernel (e.g. linear, squared exponential, Matérn, etc.) while ba = ka d') is the generic entry of the D x D coregionalization matrix B, which contains the parameters governing task-relatedness over the actions A. In the trivial case where ba,a' = we have that tasks a and a' are uncorrelated, i.c. actions and a' are unrelated the way they affect the outcome Y.
A slightly more general framework, which we adopt in this work, is given by the sum of separable ker- nels structure (Alvarez et al. 2012). This assumes that the single entry of K(-,-) reads ka,a ( v;) E9 By (x,, i.e. the sum of Q different coregionalization matrices B, with associated base kernel kg(.)). In compact matrix notation this translates into K(X,X') E Bg & Ky(X, X') for two different. input matrices X, X' E X, where & is the Kronecker product. Imposing a sum of separable kernels structure is equivalent to assuming that the collection of action-specific functions (f ())) generates from Q S D
Technically this makes the setup semi- carmmiaraaettic, rather than fully nonparametric.
4

---

Published in Transactions on Machine Learning Research (11/2022)
common underlying independent latent GP functions {ug (0)())) parametrized by their base kernel (s((s)); that is cov(u,(2;), My (2x)) = ky( 25)) (Alarrez et al. 2012
In lerms of the form of the coregionalization matrices Bg, with q E (1....,() Q), we will follow the linear model of coregionalization (LMC), which assumes that each Bg is equal B2 = Lggig with single entries Rg of 200 Rg represents the number of GP samples obtained from the same latent GP
(a,a") = EE function 9, (.).). Thus, adopting the LMC for causal learning is equivalent to assuming that correlation in the (-)a action-specific functions, modelled through the multitask kernel K(,,),, arises from Q different samples of Rg GP functions with the same kernel ());),, drawn from Q S D different independent. latent GP processes ())()) To express this more compactly, a causal multitask GP model under the K(.,))), LMC reads (-) ~ GP(m()). (((m))..((;))))) with single entries of K(:.-) being
Q AgAl (2),) implying
ka (21, 23y
Σ 9=1
By (21,
Σ
q=1
(5)
Q ΣΣ
cov(fa(ma),
cov (or(u,(ax)) of (m;)).
9.9"
In our specific case, as in Alaa & van der Schaar (2018] we employ a special case of LMC, named intrinsic coregionalization model (ICM) (Bonilla et al. 2008), where the underlying latent GP function is unique (Q = 1), so that (2,, xz) = B K (x,, x²),, with unique base kernel K(,-). The ICM specification attempts to avoid severe parameter proliferation in high-dimensional settings with multiple correlated actions D = [A], while still being capable of capturing task-relatedness through the relatively simple structure of B. However, beside the issue of parameter proliferation when A features multiple discrete actions, exact GP regression is also known to scale poorly with sample size and cardinality of input space (,8, and direct likelihood maximization methods face issues in over- over-varannn parametrized models, although some solutions, such as variational methods Titsias 2009 Hensman et al. 2013), might be adopted for better scalability.
3.2 Why multitask counterfactual learning?
We know that asymptotically, under no sample selection bias, the best approach to estimate the causal quantities fy (a:)) would be a "T-Learner" (Künzel et al., 2017 Caron et al. 2022a), which implies splitting the sample and fitting separate models for each arm A = a, f(x)a. However, in finite samples flawed by selection bias, which is the typical case of observational studies, this is often not the best strategy (Hahn) let al., 2020: Caron et al. 2022a). By simply extending the result in Alaa & van der Schaar (2018), we hereby show how increasing action and covariates spaces make the problem of estimating causal effects harder in terms of minimax risk, i.c. in terms of optimal convergence rates for any nonparametric estimator for a given space of functions. We particularly consider the problem of estimating Conditional Average Treatment Effects (CATE) between action and 6, defined under the SCM specified by (2) as Ta.)( [[Y; = a), X, a] - E[Y, (do(A, = b), X, = x = (x,) fox))) We assume that all the fa. Va E AA belong to the Hölder ball class of -mooth functions H(aa), with ag - 1 bounded derivatives and bounded in sup-norm by a constant C > 0. Considoring an no2m loss function on Ta.s (I;), namely E ((]laa,((x;) - (=.s( the difficulty of CATE estimation with a nonparametric model vi E V can be specified by the optimal minimax rate of convergence, that we define as follows (proof in the Appendix A section).
Corollary 3.1 (Minimax rate). Assume covariale space is X = [0, and Sa depends on Pa covariales s.t. P. < P, Ya A. Define na, < N as the subsample identified by action and b. If both fo E H(a.) and fo E Hl(a), then CATE optimal minimar rate for a -nor loss function is:
(1+ (1+2(Be (+= 2²))
pPoto
inf sup ([Ilta,((xs) - (Z))()) x nas Fe JuiJ's
Vlog
P pft
where x V max y) and x is asymptotical equivalence. The first terms on the RHS of relates to the problem of CATE function approximation, while the second term to the degree of sparsity of the CATE. Thus,
5

---

Published in Transactions on Machine Learning Research (11/2022)
the optimal minimax rate, which minimizes the loss in the worst case scenario permitted, is asymptotically as complex as the hardest of these two tasks. The "tightness" of this rate depends on the cardinality of the predictor space P and of the subsample defined by action a and b. This implies that, in the presence of multiple discrete actions (where some arms are likely to be very imbalancod), the rate will inevitably grows larger. For these reasons, a multitask GP prior over the actions is well suited to tackle selection bias and particularly estimation in regions with poor overlap, i.c. regions in X where we mainly observe data points with specific action A. = a and very few others. In addition to this, as shown by AAaa & van der Schaar (0118), multitask GP prior can achieve the optimal minimax rate of Corollary 3.1 in its posterior consistency.
Hence, splitting the sample into na subgroups and litting independent models can be very sample inefficient in these settings ( Multitask GPs can aid extrapolation in such cases of strong sample selection bias, by learning the correlated functions (fy, (+) jointly as fr).). The first row plots Figure 2 provide a very simple one-covariate example of how multitask lcarning addresses the issuc of extrapolation and prediction in poor overlap regions. Filting the two surfaces (3r) and (xo())) (dashed lines) through separate GP regressions results in bad fit out of overlap regions (top-left plot) in this specific case. Multitask coregionalized GP attempts to fix this problem by embedding the assumption that the two surfaces share similar patterns via joint estimation of fv. (z;) and their task-relatedness parameters, increasing sample efficiency (right panel). When the two surfaces share minor patterns instead, such as in the second row plots of Figure [ the sample efficiency gains are less siguificant; and in some more extreme cases where the surfaces do not share any similar pattern at all, assuming a multitask GP prior might also introduce bias. The issue of partial overlap might be less severe in scenarios with larger sample size, but not always; however, in settings with strong sample selection bias, scttings with multiple discrete actions or action spaces that grow with the sample size, the issue remains relevant. This is because the sampling mechanism is inherently flawed, so that even with infinite samples, poor overlap regions do not fade away, independently of the modelling assumptions one makes.
Modelling assumptions on how to estimate {sY, ()) hence ultimately depends on domain expert knowledge about the setting Hahn et al. (2020): Caron et al. (2022c). If one possesses prior knowledge that the causal effects might be fairly homogencous across units with different covariate realizations X; = 24, so that the CATE function is likely to display simple patterns (first row, Figure 2), then using multitask approach would actually help incorporate this assumption in the model. Conversely, if one believes instead that CATE is likely to be a rather complex function itself (second row, Figure 2 estimating {fr. (-)) independently would be a better choice.
3.3 Multiple Output Designs
Reverting back to setups with multiple correlated outcomes, we introduce a simple extension to the class of counterfactual GPs presented above that involves an extra multitask learning layer over the M outcomes Y,, in addition the one over AG. The way we formulate this is through an additional corogionalization matrix, in the same fashion as we did for the actions case. This exlended version featuring slacked coregionalization, has GP prior of the following form:
Y, fv. +E,, E(E,) =0 fv.0)~9P(0.Kv Ky(,)= = By & BA BABK(( k(.))
where By is the M M coregionalization matrix over the outcomes, BA the D x D coregionalization matrix over the actions and (...)) is the base kernel. The valued function fr. (-) in this case includes all the single-valued functionals (fa.m())) (.)) The extra multitask lcarning layer over the outcomes Y, is aimed at increasing sample efficiency by borrowing information among correlated outcomes, as opposed to fitting M separate counterfactual GPs with a single coregionalization layer over A, but is also conceptually sound, as the quantity of interest is indoed the joint interventional distribution p(Y|do( [do(A A = a). X = x)), which accounts for and explicitly models correlation between the outcomes, rather than the collection of marginal distributions {p(Y... |do(A A = a)) m= which leaves correlaiin unspecified. As we will address in the later
6

---

Published in Transactions on Machine Learning Research (11/2022)
Rt
Ft
Treeted Cuntrol
Treated Control
Figure Simple one covariate example, with {0, 13. Overlap is guarantcod to hold over the whole support X' in the data generating process, i.e. every unit has non-zero probability of being assigned to either A, =1 A, = O, but p(A, 1[x,) is generated as increasing function of X, (selection bias). In the top row simulation, the two underlying counterfactual surfaces fYa (z44 (dashed lines) are generated with very similar patterns, thus GP (left panel) is unable to borrow information from the other arm in poor overlap regions contrary to multitask GP (right panel). In the bottom row simulation instead we generate less similar surfaces, so borrowing of information through multitask GP does not lead to any improvement.
section, although the extra layer defined by By allows for higher sample efficiency, it also poses some issues due parameter proliferation and stability of the optimization problem in high dimensions.
4 Counterfactual Multitask Deep Kernel Learning
Gaussian Processes regressions known to scale poorly with high dimensions. Their typical computational cost amounts to (n²) ) for training points and O(n²) for test points. Similarly, coregionalized GPs suller from prer-parametrization and instability in the optimization procedure as the mumber of inputs P and the number of discrete actions increase. Deep Kernel Learning (DKL) was firstly introduced by Wilson et al. (206)) with the aim of combining the scalability of Deep Learning methods and the nouparametrie Bayesian uncertainty quantification of GPs tackling prediction tasks in high-dimenssonal settings. Given a base kernel k(:--) (e.g, linear, squared exponential, etc.), DKL structure learns functional fy. (-) by passing the P inputs X; X' through deep architecture (a colly-coneeete feedforward neural network in our case), which maps them to lower dimensional representation space via non-linear activation functions. The base kernel k(---) is then applied in this lower dimensional representation space, ((h(()) ((h)))))))))) where (()) is a neural network's s hidden layer, constituting final Gaussian Process layer (or an infinite basis functions representa a tion layer). The resulting mathematical object can be described as a kernel being applied to a concatenation of linear and non-linear functions of the inputs, namely K(x,, x²) = K(g ....... o gi (2:), 91 G ... o (xx)) (Bohn) et al. 2019). Thus, the DKL architecture is end-to-end, fully-connected and learnt jointly: the P inputs are passed to l hidden neural nets layers where the last hidden layer before the GP layer typically maps them to lower dimensional representation space (with e.g. two hidden units). This is what generates benefits in terms of scalability compared to classic GP, as the base kernel k(.,-) is applied to a lower dimensional representation space, rather than the higher dimensional inputs space directly. Another intrinsic ad advanlage of DKL is that the deep architecture preceeding the GP layer can itself learn arbitrarily complex function, so the choice of specific GP kernel becomes less cunbersome. For example, Wilssnn et al. (2016) show that DKL is more robust in recovering step functions, due to weaker smoothness assumptions compared standard GP kernels.
7

---

Published in Transactions on Machine Learning Research (11/2022)
x Japer
Tegect leyer
Hibden byarss
Output beyer
A(4) ( b1).
A() s, 611- + by).
f(h1)) + [6] gp(O, R(-.)))
R(..)= @ BA Ko(h
Figure 3: Counterfactual multitask DKL architecture. The P raw inputs are passed through a deep learning structure with hidden layers. Multioutput separable kernels (inducing coregionalization over actions A and outcomes Y) are then applied to the last Gaussian Process hidden layer, before the M action-specific output layer. Parameters are estimated jointly by minimizing the negative log likelihood.
DKL naturally presents some limitations concerning the more burdensome parameter tuning (e.g., hidden layers and units selection) and the fact that they more easily tend to overfit when avely-paramettixxe (we refer to Ober et al. (2021) for a more detailed discussion of the issue). The kernel .) in the last GP layer of DKL architecture can easily incorporate the separable kernel structure for multitask learning, in the same fashion as the class of causal GPs presented earlier. Thus, we propose multitask modelling framework to induce correlation across the action-specifie functions {fy, (-)) under the name of Counterfactual DKL (CounterDKL)), where similar Intrinsic Coregionalization Model (ICM) in the same spirit of ( is placed the last hidden layer of the neural network (0) 9: (t-1) 4 bi), where (Wr, bı) are the last layer's weights and bias, such that
fy. (h gP (P(0.((..)) (Th",,0' KA(a, a') & k(h") ,00) B (10,,,000)
where
(8)
In this case, the Kronecker product of the coregionalization matrix occurs the last hidden layer, and features lower dimensional representations instead of the potentially large number of raw inputs. Similarly, can induce coregionalization over the M outcomes by adding another level of corogionalization, with the kernel reading ((h"), k0') (u, & KA(a, a') & K,(,",, h(") By & BA @ K, (1,9), n.0).. Figure graphically depicts counterfactual multitask Deep Learning architecture, with fully-connected hidden layers, final (infinite) GP layer and the M action-specific outcomes. The parameter space in CounterDKL comprises: i) the set of deep neural network's weight matrices (W.).-1 and biases {b. ii) the base kernel K hyperparameters $, e.g. in the case of squared exponential kernel & includes lengthscales and variances parameters & (e, iii) the coregionalization matrix B entries. Hence, the parameter space is the collection θ ((W,) {b, \/ Φ B). These parameters are estimated jointly via maximization of the log- oog-marginal likelihood. More details are provided in the Appendix Section D and in Wilson et al. (2016): Gardner et al. [2018). In the next section, we will investigate properties of counterfactual multitask GPs and DKL on a variety of experiments.
5 Experiments
The fundamental problem of causal inference is that the interventional quantity (d([ = a). X; = x:)) is never observable, so we have to resort to simulation to fully evaluate the methods on individual causal effects (ICE) estimation. We evaluate the performance of counterfactual GPs and counterfactual DKL on data generating process with three different tasks, and on a worldworld example combining experimental and observational data. For the first simulated experiment, construct the DGP such that the backdoor
8

---

Published in Transactions on Machine Learning Research (11/2022)
Varying F-10 OPE
Gr
2.0
3.
0.30
DaL
Moot
MEDE.
0.73
12-
0.70
10-
510
3000
1500
2000
2003
570
2700
1500
2000
2500
530
2090
3500
2002
2500
Verying No.1500
0.60
0.32
0.13
14.
0.70
12
10-
08
15
Figure 4: Results on performance of the methods compared, in terms of RMSE or Optimal Allocation Rate (OAR), averaged across B 100 replications for each N E {500, 1000, 1500, 2000, 2500} (first row) and each P E (10, 15, 20, 25} (second row). First column: RMSE evaluated on the individual causal effect (TCE) estimation task (on the test set). Second column: RMSE evaluated on the OPE task. Third column: OAR on the OPL lask, defined as percentage of units correctly allocated to the best action among the D ones.
crilerion holds for X, E X' The GPs and DKL implementations in the simulated examples all make use of the KISS-GP approximation to compute the base kernel covariance matrix as K2 = M MK M in the GP 0,0
layer for better scalability (Wilson & Nickisch 2015 Wilson et al. 2016 4
5.1 Simulated Example
We consider simulated setting with D = 4 possible actions A= (0, and M = 2 correlated outcomes Y = (Y1,Y)) E RR.. Actions and outcomes are generated according to a policy Tb( = p(A, = a|X; = x;) and an outcome function (x,), both dependent on the covariates X; E X. The probabilistic DGP is fully described in the Appendix B of supplementary materials. The models we compare are the following: i) separate standard GP regressions, employed to fit (.,()) distinctly for each outcome and for each action (GP); ii) counterfactual multitask GP regression (Alaa & van der Schaar) 2017 2018), with coregionalization over A only, meaning that fit two separate models for each outcome, but unique multi-valued function model for A, (CounterGP)) ); ili) counterfactual multioutput GP regression, a unique model with coregionalization both over A; and Y; (MOGP): iv) separate DKL regressions with 3 hidden layers of [50, 50, 2] units, the equivalent of i) but with deep kernel implementation (DKL); v) counterfactual multitask DKL regression with 3 hidden layers of 50, 2] units, the DKL equivalent of ii) (CounterDKL); vi) counterfactual multioutput DKL, the DKL equivalent of îii) (MODKL). In particular, we consider two slightly different versions of this setup. In the first version we fix the number of covariates to p= 10 (only 7 of them being relevant for the estimation) and study the behaviour of the estimators with increasing sample size N (500, 1000, 1500, 2500). Tn the second version we fix sample size to N = 1500 and study the behaviour of the estimators with increasing number of covariates P E {10, 15, 20, 25}. Performance of the models is evaluated on the following three related tasks:
ICE: The first is the prediction of Individual Causal Effects (ICE). This tackles the estimation of the average causal effect of playing action A, outcome Y,, given certain realization of the covariates
4 All experiments were run on Intel(R) Core(TM) 17-7500U CPU @ 2.70GHz, 8Gb RAM CPU.
9

---

.5, 25). Performance is Published in Transactions on Machine Learning Research (11/2022)
1.0
1.00
CuaranaP
1.6.
0.55
MOBCI
1.4
Coverage
1952
1.2
0.90
Y:Degeee Contounding
Degree Contounding
Figure 5: Models' performance in terms of RMSE (left plot) and 95% Coverage (right plot), in estimating Individual Causal Effects (ICE) on a 20% left-out test set, given an increasing level of confounding, representod by the Y parameter: higher values of y corresponds to higher probability of being assigned to of the two action A; f3, 4), thus generating more arms imbalance.
space, X; = x, i.e. the estimation of ICE: B(Y;)do(A, - a) X: wx). This is carried out using a 80% training set, and evaluated via RMSE on 20% left-out test set.
OPE:
The second is Off-Policy Evaluation, which is concerned with quantifying how good a given alternative policy ne is, compared to the action allocation policy that generated the data (bchavior policy, TH). This done by estimating the policy value, defined as the cumulative reward V(TT) = a:)]] originating from me. In our case we pick the alternative policy action allocation, i.e. Te Multinom(2255
E Ex,A. EE Te (a.l)x)) (Yildo(A)
to be a uniformly-at-randee evaluated through RMSE the entire sample.
The last is Off-Policy Learning, which is concerned with finding the optimal policy x² : + A, that is the policy that generates the highest policy value: qे. E arg max V(Tp)) This last task is evaluated through an Accuracy metric on the whole sample, which we label Optimal Allocation Rate (OAR), indicating the percentage of units correctly assigned to their specific optimal action H² (xx) = a, i.c. the action that generates the best outcome for them.
OPL:
Since we are dealing with M = outcomes, we produce performance measurements on RMSE and optimal allocation rate for both outcomes and then average them, assuming both outcomes are given equal policy importance and live on the same scale. For both versions of the setup, namely increasing N and increasing P, replicate the experiment B 100 times to obtain Monte Carlo averages and 95% confidence intervals for the metrics. Results are depicted in Figure g RMSE performance in all models for increasing N and P behaves accordingly to Corollary 3.1 CounterDKL. and MODKL. perform consistently better than the GP models, as they scale better with an increasing sample size N and increasing number of predictors P. Particularly MOGP's performance deteriorates for issues related to stability of the marginal likelihood maximization and over-parameetriztiion, as we had to omit. it from the study of increasing predictors due to failed convergence for P 10. The advantages over standard DKI. regression instcad are entirely attributable to sample efficiency gains from multilask coregionalization in CounterDKL and MODKL, both in the increasing N and increasing P studies. In the case of increasing P, emphasize that as predictor space grows larger the causal DGP becomes relatively sparser (only 7 predictors out of P remain relevant for the estimation), especially in the case of P = {20, 25). So in these two cases the batch of DKL models would perhaps achieve better performance from increasing the number of hidden units or hidden layers and adding regularization (dropout, εl or /2 regularizers) in the deep architecture part.
In addition, we run slightly different version of the ICE experiment above, to further investigate properties of the models in terms of uncertainty quantification, that we measure through the 95% coverage of each ICE
10

---

Published in Transactions on Machine Learning Research (11/2022)
Model
Train MAE
Test MAE
Train Rpol
Test Rpol
Runtime (s)
GP CounterGP
0.033 + 0.006
0.036 + 0,008
0.22 0.02
0.27 + 0.02
171.3 + 16,1
0.033 + 0.006
0.035 + 0.007
0.24 0.01
0.27 + 0.02
248.6 + 6.4
PCA + GP
0.073 = 0.002
0.074 = 0.003
0.22 0.01
0.27 = 0.02
66.3 + 2.4
PCA + CounterGP
0.074 + 0.001
0.074 +: 0.001
0.23 +: 0.01
0.26 + 0.02
126.1 + 3.9
AutoEnc + GP
0.075 = 0.004
0.075 = 0.003
0.21 0.03
0.27 = 0.02
76.0 + 3.0
AutoEnc + CounterGP
0.076 = 0.003
0.076 = 0.003
0.24 = 0.02
0.30 + 0.03
138.7 + 9.2
DKL
0.029 + 0.011
0.042 + 0.015
0.20 + 0.01
0.21 + 0.02
44.8 + 3.3
CounterDKI.
0.011 = 0.003
0.015 + 0.005
0.22 0.01
0.25 + 0.02
122.7 + 7.4
Table 1: Train and test set performance on the Jobs data experiment in terms of Mean Absolute Error (MAE) in estimating ATT, Policy Risk (Roo)) and overall runtime (s), with 10-fold cross-validated 95% intervals. Bold indicates best performance.
estimates (then averaged over actions and outcomes). This is defined as
Coveragegs% = N Σ 1 (fa(z)))on) S fa(x;) S fa(wo))mo))
where (2))))o and fa(()) Jupp are the lower and upper bands of the 95% credible interval output by the model on fo (2,), while fxx) is the true individual counterfactual outcome correspouding to action A, = a. Given fixed N 2000 and parameter 7, which governs the level of confounding, or the degree of groups imbalance in terms of action allocation. Particularly, for increasing values of y, we assign higher probability of choosing one of the two action A, = {3, 4). This generates action arms imbalance as it leaves gradually less units in arms { 2}. Results are gathered in Figure 5 where MODKL display higher performance both in terms of error and uncertainty quantification.
20 and similar data generating process as before, we introduce the
5.2 Real-World Example: Job Training Programs and Unemployment
We demonstrate the efficiency of CounterDKL, also on a second experiment taken from Shalit et al. (201)) involving a popular real-world study on a job training program, dating back to LaLonde (1986). The distinctive feature of this dataset is that it combines randomized and an observational subgroups, where the aim is estimate the effects of participation on job training program on earnings and employment. The randomized experiment features 297 treated and 425 control units; The observational subsample is instead made of 2490 control units only. The binary treatment A; E {0,1} denotes participation to the job training program. The original outcome YÚ is earnings after the program, which is censored continuous (Y; = for unemployed units). However, following Shalit et al. (2017). we construct a binary indicator Y E 10, denoting employment status at the end of the job training program as outcome. This gives us the opportunity to demonstrate the use of the methods presented in this paper also on binary/categore type of outcomes. To this end we use the classification method for GPs proposed in Milios et al. (2018), where class labels are interpreted as coefficients of degenerate Dirichlet distribution, which makes the GP classification task efficiently faster and more scalable. The 7 covariates X, E X in the study are the following: ago, years of schooling, african american ethnicity, hispanic ethnicity, marital status, high school diploma. Given the presence of a randomized subsample, we can exploit it to compute unbiased estimates of the quantities of interest and treat them as ground truth. The two quantities of interest in this case are: i) the Average Treatment Effect on the Treated group (ATT), defined as ATT = 7- Σle 9: C-1 E Vis where T and C are the number of Ireated and control units in the experimental data; ii) the Policy Risk (Shalit et al.) 2017), defined as the average error in allocating the treatment according to the ICE estimates policy rule namely <xx;)) = if ICE (Y(] = 1), 2,) do(,, = 0), m;) Rpol- (EYYddo(AA 1), TZx,) 1)p(=(<z,) 1) + E(Y|do(A 0), TT.,) = ()p(=(z;) -)]. Notice that cannot measure performance ICE directly this is always unobservable in real-world scenarios;
or
11

---

Published in Transactions on Machine Learning Research (11/2022)
also, restrict analysis of average causal/treatment effects the treated group since are sure that overlap holds there, as all the treated units were part of the randomized experiment subgroup, while the observational subgroup is made only of control units. More details about this experiment can be found in the Appendix C of supplementary materials and in Shalit et al. (0117). We compare the following models: i) GP and CounterGP, as in Alaa & van der Schaar (2017); ii) vanilla PCA plus either GP or CounterGP; iii) vanilla deep AutoEncoder plus either GP or CounterGP; iv) DKL and CounterDKL (ours). Results on performance are gathered in Table [ in terms of 70%-30% train and test set Mean Absolute Error (MAE) on ATT, Policy Risk Reol and average runtime, accompanied by 10-fold cross-validated 95% error intervals. In this example multitasking is induced only over the binary treatment, As the results depict, by operating jointly via a unique loss function, CoLnterDKL is significantly more efficient than naively applying dimensionality reduction and fitting multitask GP on a lower dimensional space as two separate steps. It also displays gains over CounterGP, thanks to its deep component that guarantoes better computational time (in terms of runtime) and scalability, and able learn arbitrarily complex functions while imposing weaker smoothness assumptions than standard GP kernels, even on a low-dimensional covariate space example such as the one presented here (7 covariates).
we deal with just a single outcome of interest.
5.3 The Infant Health Development Program data
Finally we compare CounterDKL with few other recent methods for causal effects estimation, on popular simulated experiment utilizing the Infant Health De-
Model
Train RMSE
Test RMSE
velopment. Program (THDP) data, originally found in
RF
1.85 + 0.13
2.39 t 0.17
Hill 2001), and more recently in several contributions on Conditional Average Treatment Effects (CATE) es-
X-RF
3.29 1 0.23
3.37 1 0.24
CF
3.10 + 0.21
3.07 1 0.20
timation, such as Shalit et al. (2017); Alaa & van der Schaar (2017 2018); Caron et al. (2022a) among oth- ers. The experiment is a semi-simulated setup. It makes use of real-world data from the IHDP study, randomized clinical trial aimed at improving the health status of premature infants with low weight at birth through pediatric follow-ups and parent support groups, and recreates an observational type of study by removing a non-random portion of treated units, namely those with "non-white mothers". This leaves total of 139 observations in the treated group and 608 in the control. In addition, the semi-simulated setup uses the real-world binary treatment A, {0, 1} and the 25 available covariates X; E X, but simulates the two continuous potential outcomes (Y1 Yo) E R as described in the non-linear "Response Surface B" setting in Hill (201)). As anticipated above, the estimand of interest in this case is the Conditional Average
BART
0.97 0.04
1.44 1 0.09
X-BART
2.07 + 0,14
2.13 0.14
BCF
0.87 0.05
1.34 0.10
CounterGP
0.61 0.02
0.70 + 0.04
CounterDKL
0.62 + 0.02
0.67 + 0.04
Table 2: RMSE of compared models on the semi- simulated THDP setup, evaluated for CATE estimation on 80%-20% train-test sets.
Treatment Effects (CATE), defined as |[ [d([, = 1), = x) E [do( = 0), X; = xx). The models we compare include: i) vanilla Random Forest (RF), as a T-Learner (K Künzel et al. 2017 Caron et al., 2022a); ii)
X-Learner version of Random Forests (X-RF) as in (Künzel et al. (017): iii) Causal Forest (CF), or Random Forests as an R-Learner, developed by Wager & Athey| (2018); iv) vanilla Bayesian Additive Regression Trees (BART), as a T-Learner; v) X-Learner version of BART (X-BART); vi) Bayesian Causal Forests (BCF)
by Hahn et al. (2020); vii) Counterfactual GP (CounterGP) as in Alaa & van der Schaar (0018): ; viii) our
Counterfactual DKL (CounterDKL) with [100, 100, 2] hidden layers. Results reported in Table refers to 1000 replication of the experiment on 80%-20% train-test split as in Alaa & van der Schaar 2017
6 Conclusions
Throughout this work, we considered the problem of counterfactual effects learning using observational data, which is of interest in domains where exploration of policies is costly (healthcare, socio-economic sciences, etc.). We reviewed the class of counterfactual GP regression models, extending to adjust to multiple actions and outcomes settings, and discussed how multitask learning over multiple actions helps addressing finite sample selection bias. We then introduced new class of counterfactual models based on
12

---

Published in Transactions on Machine Learning Research (11/2022)
Deep Kernel Learning, whose main advantages lie in their more flexible function approximation capabilities and better scalability. While counterfactual GPs struggle to scale up with sample size, number of predictors and number of actions/outconm to coregionalize over, DKL capitalizes on these components by learning lower dimensional representations. We stress that the class of DKL methods proposed can be easily expanded to carry out counterfactual learning in other more complex scenarios such as: i) unobserved confounding where identification is still possible (instrumental variables); ii) dynamic settings such as dynamic treatment regimes or RL; iii) non-standard data like images (as DKL can incorporate any type of deep architecture).
References
Ahmed Alaa and Mihacla van der Schaar. Limits of estimating heterogencons treatment effects: Guidelines for practical algorithm design. In Proceedings of the 35th International Conference on Machine Learning, pp. 129 138, 2018.
Ahmed M. Alaa and Mihaela van der Schaar. Bayesian inference of individualized treatment effects using multi-task Gaussian Processes. In Proceedings of the 3ist International Conference on Newral Information Processing Systems, NIPS'17, pp. 1427-, 2017.
Joshua D. Angrist, Guido W. Imbens, and Donald B. Rubin. Identification of causal effects using instrumental variables. Journal of the American Statistical Association, 91 91443)):444455, 1996.
Susan Athey and Stefan Wager. Policy Learning With Observational Data. Econometrica, 89(1):133 161, January 2021.
Heejung Bang and James M. Robins. Doubly robust estimation in missing data and causal inference models. Biomelrics, (4):9962 973, 2005.
Elias Bareinboim, Andrew Forney, and Judea Pearl. Bandits with unobserved confounders: A causal approach. In Advances in Newral Information Processing Systems 29, volume 28, 2015.
Bastian Bohn, Christian Rieger, and Michael Griebel. A representer theorem for deep kernel learning. J. Mach. Learn. Res., 20(1):2302 2333, 2019.
Edwin V Bonilla, Kian Chai, and Christopher Williams. Multi-task Gaussian Process prediction. In Advances in Neural Information Processing Systems, volume 20, 2008.
Alberto Caron, Gianluca Baio, and Ioanna Manolopoulou. Estimating individual treatment effects using non-parametric regression models: review. Journal of the Royal Statistical Society: Series A (Statistics in Society), 2022a. doi: https://abi.orr//..11111 12824.
Alberto Caron, Gianluca Baio, and Ioanna Manolopoulou. Shrinkage Bayesian Causal Forests for heterogeneous treatment offects ostimation. Journal of Computational and Graphical Statistics, pp. 1-13, 2022b. doi: 10.1050/10618600.2000 20222467999....
Alberto Caron, Gianluca Baio, and Toanna Manolopoulou. Interpretable deep causal learning for moderation effects. In ICML 2022, 2nd Interpretable Machine Learning for Healthcare Workshop, 2022c. URL
https:///rxiv. org/ waws//200 10261
Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins. Double Double/debiasee debiased machine learning for treatment and structural parameters. The Econometrics Journal, [1(1)C11-666 2018.
Miroslav Dudík, John Langford, and Lihong Li. Doubly robust policy evaluation and learning. In Proceedings of the 28th International Conference on International Conference on Machine Learning, pp. 1097 1104, 2011.
Miroslav Dudik, Dumitru Erhan, John Langford, and Lihong Li. Doubly robust policy evaluation and optimization. Statistical Science, 29(4):485-51 2014.
13

---

Published in Transactions on Machine Learning Research (11/2022)
Mehrdad Farajtabar, Yinlam Chow, and Mohammad Ghavamzadch. More robust doubly robust off-policy evaluation. In Proceedings of the 35th International Conference on Machine Learning, volume 80, pp- 1447 1456, 2018.
Jacob R. Gardner, Geoff Pleiss, David Bindel, Kilian Q. Weinberger, and Andrew Gordon Wilson. Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration. In Proceedings of the 32nd International Conference on Newral Information Processing Systems, pp. 7587 7597, 2018.
P. Richard Hahn, Jared S. Murray, and Carlos M. Carvalho. Bayesian Regression Tree Models for Causal Inference: Regularivation, Confounding, and Heterogeneous Effects. Bayesian Analysis, 15(3):965 1056, 2020.
Jason Hartford, Greg Lewis, Kevin Beyton-Browar and Matt Taddy, Deep IV: A flexible approach for counterfactual prediction. In Proceedings of the 34th International Conference on Machine Learning, volume 70, pp. 1414-1423, 2017.
James Hensman, Nicolò Fusi, and Neil D. Lawrence. Gaussian processes for big data. In Proceedings of the Twenty-Ninth Conference on Uncertoinly in Artificial Intelligence, UAT13, pp. 282-290, 2013.
Jennifer L. Hill. Bayesian nonparametric modeling for causal inference. Journal of Computational and Graphical Statistics, 20(1):217 240, 2011.
R. Hodson. Precision medicine. Nature, 577(9),,), 2016.
D. G. Horvitz and D. J. Thompson. A generalization of sampling without replacement from finite universe. Journal of the American Statistical Association, 47(260):663-66) 685, 1952.
Kosuke Imai and David A van Dyk. Causal inference with general treatment regimes. Journal of the American Statistical Association, 99(467):854 866, 2004.
Guido W. Imbens and Donald B. Rubin. Causol Inference for Statistics, Social, and Biomedical Sciences; An Introduction. Cambridge University Press, 2015.
Nan Jiang and Lihong Li. Doubly robust off-policy value evaluation for reinforcement learning. In Proceedings of The 33rd International Conference on Machine Leorning, volume 48, pp. 652-661, 2016.
Fredrik Johansson, Uri Shalit, and David Sontag. Learning representations for counterfactual inference. In Proceedings of The 33rd International Conference on Machine Learning, volume 48, pp. 3020-3029, 2016.
Nathan Kallus. Balanced policy evaluation and learning. In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pp. 8909 8920, 2018.
Nathan Kallus. More efficient policy learning via optimal retargeting. Journal of the American Statistical Association, 116(534):646--4 (68866666666666 2021.
Toru Kitagawa and Aleksey Tetenov. Who should be treated? empirical welfare maximization methods for treatment choice. Econometrica, 86(2):591-66 616, 2018.
Sören Kiinzel, Jasject Sekhon, Peter Bickel, and Bin Yu. Meta-learners for estimating heterogencous treatment effects using machine learning. Proceedings of the National Acndemy of Sciences, 116, 06 2017.
Robert J. LaLonde. Evaluating the econometrie evaluations of training programs with experimental data. The American Economic Review, 76:604 620, 1986.
Dimitrios Milios, Raffaello Camoriano, Pietro Michiardi, Lorenzo Rosasco, and Maurizio Filippone. Dirichlet- based gaussian processes for large-scale calibrated classification. In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pp. 6008-6018, 2018.
Whitney K. Newey and James Powell. Instrumental variable estimation of nonparametric models. Econometrica, 5165::1565-1555 -178, 2003.
14

---

Published in Transactions on Machine Learning Research (11/2022)
X Nie and S Wager. Quasi-oracle estimation of heterogeneous treatment effects. Biometrika, 108(2):299-33)) 09 2020.
Xinkun Nie, Emma Brunskill, and Stefan Wager. Learning when-to-treat policies. Journal of the American Statistical Association, O(ja):1 58, 2020.
Sebastian W. Ober, Carl E. Rasmussen, and Mark van der Wilk. The promises and pitfalls of deep kernel learning. In Proceedings of the Shirth-Seventh Conference Uncertainty in Artificial Intelligence, volume 161, pp. 1206-1216, 2021.
Judea Pearl. Causality: Models, Reasoning and Inference. Cambridge University Press, USA, 2nd edition, 2009. ISBN 052189560X.
Min Qian and Susan A. Murphy, Performance guarantees for individualized treatment rules. Ann. Statist., 19(2):1180-1210, 04 2011.
Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian Processes for Machine Learning (Adaptive Compulation and Machine Learning). The MIT Press, 2005.
Paul R. Rosenbaum and Donald B. Rubin. The central role of the propensity score observational studies for causal effects. Biometrika, 70(1):41 55, 04 1983.
Donald B. Rubin. Bayesian inference for causal effects: The role of randomization. Ann. Statist., 6(1):34 58, 01 1978.
Phillip J. Schulte, Anastasios A. Tsiatis, Eric B. Laber, and Marie Davidian. Q-and a-learning methods for estimating optimal dynamic treatment regimes. Statistical Science, 29(4):640-66 661,666,111 2014.
Cri Shalit, Fredrik D. Johansson, and David Sontag. Estimating individual treatment effect: Generalization bounds and algorithms. In Proceedings of the 34th International Conference Machine Learning Volume 70, volune 70, pp. 3076 3085, 2017.
Richard S. Sutton and Andrew G. Barto. Reinforcement Learning: An Introduction. A Bradford Book, 2018.
Yee Whye Teh, Matthias Seeger, and Michael I. Jordan. Semiparametric latent factor models. In Proceedings of the 10th International Workshop on Artificial Intelligence and Statistics, volume R5, pp. 333 533-340, 2005.
Philip Thomas and Emma Brunskill. Data-efficient off-policy policy evaluation for reinforcement learning. In Proceedings of The 33rd International Conference on Machine Learning, volume 48, pp. 2139-2148, 2016.
Jin Tian and Judca Pearl. A general identification condition for causal effects. In Eighteenth National Conference on Artificial Intelligence pp. 567 573, 2002.
Michalis Titsias. Variational learning of inducing variables in sparse gaussian processes. In Proceedings of the 12th International Conference on Artificial Intelligence and Statistics, volume 5, pp. 567-574 2009.
Masatoshi Uehara, Masahiro Kato, and Shota Yasui. Off Off-po evaluation and learning for external validity under a covariate shift. In Advances in Neural Information Processing Systems 33, 2020.
Stefan Wager and Susan Athey. Estimation and inference of hetcrogencous treatment effects using random forests. Journal of the American Slatistical Associalion, 113((23):1225 1242, 2018.
Andrew Gordon Wilson and Hannes Nickisch. Kernel interpolation for scalable structured gaussian processes (P). In Proceedings of the 32nd International Conference on International Conference on Machine Learning, pp. 1775 1784, 2015.
Andrew Gordon Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric P. Xing. Deep kernel learning. In Arthur Gretton and Christian C. Robert (eds.), Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, volume 51, pp. 370-378. PMIRR 2016.
15

---

Published in Transactions on Machine Learning Research (11/2022)
Liuyi Yao, Sheng Li, Yaliang Li, Mengdi Huai, Jing Gao, and Aidong Zhang. Representation learning for treatment effect estimation from observational data. In Advances in Neural Information Processing Systems 31, pp. 2633 2643, 2018.
Baqun Zhang, Anastasios A. Tsiatis, Eric B. Laber, and Marie Davidian. A robust method for estimating optimal treatment regimes. Biometrics, 68(4):1000-111; 2012.
Yingqi Zhao, Donglin Zeng, A. John Rush, and Michael Kosorok. Estimating individualized treatment rules using outcome weighted learning. Journal of the American Statistical Association, 107(499): 1118, 2012.
Xin Zhou, Nicole Mayer-Hambleu, Umer Khan, and Michael R. Kosorok. Residual weighted learning for estimating individualized treatment rules. Journal of the Americon Statistical Association, 112(517): 169-187, 2017.
Mauricio A. Álvarez, Lorenzo Rosasco, and Neil D. Lawrence. Kernels for vector-valued functions: A review. Foundations and Trends in Machine Learning, 4(3):195 266, 2012.
16

---

Published in Transactions on Machine Learning Research (11/2022)
A Proofs and Discussion
In this first section of supplementary material we provide assumptions, proofs and brief discussion of the two theoretical results in the main paper (Section 2 theorem and Section 3.2 corollary).
A.1 Backdoor Adjustment
The backdoor adjustment thcorem (Paar 2009) is well known and understood results in causal inference.
Its extension to multi-action and multi-outcome settings is trivial. We briefly reformulate the necessary assumptions as follows.
Assumption A.1 (Backdoor Criterion). Using the causal graph terminology, we denote by pa(V,) the parents of a specific random variable V, and by de(V;) its descendants. We assume that pa( (A,Y, Y) S X, de( Y) & X. Namely X contains (but does not necessarily coincide with) all the common parent variables of A E A and Y [Ymm E,, for all E {I,..., MM, and does not contain any common descendant of them.
Assumption A.2 (Overlap). Defining the propensity score as (x,))) = p(A, a|X,, æ;), we require that Na (2;) E (a, a) for all E (1,..,N) N), where a E (0, 1). The scalar a represents how tight we require the overlap between units in each arm terms of the covariates to be.
Proof of Backdoor Adjustment Theorem Under Assumption 1.1 and 1.2, we are able to identify the effect + Y = {Ym mm in the form the joint interventional distribution p( do(AA = a)), as:
-f |(Y do(A= - a), do(A = a))
P do(A a))
x
/ p(Y a. x))(( / do(A a)) dæ
X
/ p(Y | ()(())
The above derivation refers to the marginal interventional distribution, with respect to X. Often we are interest in the conditional interventional distribution p( do(A= = a),X) X) instead (e.g., conditional on patient's characteristics), such as when estimating CATE: r(x) = (x) fo())) The only additional requirement
compared to the original version of backdoor adjustment (Pearl 2009 is that Assumption 1.1 holds for all
the collection of outcomes Y.
A.2 Corollary 3.1
Corollary 3.1 is a trivial extension of the result on CATE optimal minimax rate derived in Theorem 1 by Alaa & van der Schaar (2018) discrete multi-action set domains, where specifically f0, 1} A. Thus the proof follows straightforwardly from Alaa & van der Schaar 2018 The main difference is the following.
Proof of Corollary 3.1 Optimal minimax rate Alaa & van der Schaar (2018) define the hardness of CATE estimation between binary action a, b A: (a, = fa( - (xx)). If we have multi-actions space, it means that, given the whole sample size is N, each pair (strictly more than one pair) of discrete actions a, b defines a subsample (and thus a subpopulation) of na.5 < N units. This implics that the hardnoss of approximating a function in the Hölder ball class and of performing variable selection is proportional to the smaller subsample Ra,b not N, which makes the CATE estimation problem harder the smaller na,b-
This likely happens when multi-action scenarios feature infrequently explored action arms. Thus, technically speaking, result in Theorem in Alaa & van der Schaar (0088) is a special case of Corollary 3.1 where A = (0, 1).
17

---

Published in Transactions on Machine Learning Research (11/2022)
B Data Generating Processes
We hereby describe the causal data generating processes the simulated examples of the paper (Section 3.2 and Section 5.1).
B.1 Section 3.2 one covariate example
For the simple one-covariate example in Section 3.2 (Figure 2), where discuss the benefits of multitask counterfactual learning, generated N 300 data points from one, uniformly distributed covariate, X; ~ Uniform(-3.3) Then we generated a binary action variable A, ~ Bernoull(p(4: p(A, 1|z,) ( and (--) is the standard normal edf. Finally, the two counterfactual outcome surfaces were generated as fo(=;) = 2 + 0.3 X, and fi(2;) = 3 fo(x;), with the final outcome being Y (o( + A,; + E, where T(x;)= = (144) (x:))) is the CATE function and é, ~ (0, 0.75²).
1 [x:)), where
B.2 Section 5.1 experiment
The causal data generating process for the simulated experiment of Section 5.1 described as follows. The P covariates are generated from a uniform distribution X Unif(-3,3) for j (I....P) and NJ.. The action allocation policy is simulated according to multinomial distribution where the probabilities of being assigned to action A, are generated softmax function of the covariates ΣeEA exp(X,Ba)) where 8a is an action-specific P-dimensional sparse
p(A; al al =.) vector of action-specific coefficients defined as follow:
8 = [-1
-0.8
-0.1
-011
0
0].
B2 0
0
0.8
0.2
0
" 0]
B3 115
-0.8
0.1
0.1
0 ... 0]
B. [-1
-0.8
0.1
0.1
0]
Thus A, drawn from a multinomial with vector probabilities parameter p(A, = = =x)) The A 2 action-specific correlated counterfactual outcomes Y, do(A; = a) instead are generated as
Y do(A; a) fr. + tes E; ~ N(O, E..),
where:
frn = 3 + 0.4X,, - 0.3X; + 0x2exp(Xx) + 0.6 (X,,
fyi2 11 fri 0.1X,
fY 1 ++ fru, 0.3X,
fYl 0.5 + frn 0.5X,
fy21 = 1 + 0.2X, - 0,2X3 0x1xxp(Xx)
fY22 = 2 fr2 + 0.2X3
SY23 2 fy2 + 0.4X,
fY24 1 fr2ı + 0.5X6
and where diag(E.) [o1,..., 04]. with 01 ... 04 0.5, and off-diagonal elements are 0. Finally, briefly describe the main specifications of the methods compared. The GP models (GP, CounterGP and MOGP) all employed a RBF base kernel, while the DKL models employed three [50, 50, 2] hidden layers feedforward neural network before the GP x-layer, which itself employs a RBF base kernel. The multitask and multioutput models (both GPs and DKLs) all make use of the Intrinsic Coregionalization Model (ICM), such that K(x,, = By Ba & (2,,,, a²).. ). All model were optimized through the Adam solver. More details and fully reprodueible code on this experiment can be found in the Github repository: https://github. http://ggtthbbb.coaaaa/mma http://////ccccco//ooouuui Coo.aterr/n/
18

---

Published in Transactions on Machine Learning Research (11/2022)
GP
CounterGI
DKL
CounterDKL
Data
N
P
# actions
indian
0.390
0.392
0.376
0.347
indian
573
heart
2.553
1.076
0.433
0.410
heart
270
13
2
yeast 0.534
0.657
1.3144
0.081
yeast
1484
10
contracep
0.339
0.003
0.008
0.007
contracept
1473
3
Table 4: OPE absolute regret on UCI datasets. Bold denotes best performance.
Table 3: UCI datasets characteristics.
C The Job Training Data
The Job Training data (aLLonde] 1986) are a popular case study in the causal inference literature. They comprise a portion data pertaining to a randomized experiment and portion of observational data, with the randomized exporiment featuring 297 treated and 425 control units, while the observational data being of 2490 control units only. Given the randomized subsample of the data, we can obtain an unbiased estimate (computed the randomized units only) for the Average Treatment Effect the Treated group (ATT) as ATT = T 1/ C-1 Ei= Vi, where T and C' are the number of treated and control units in the experimental data, and treat this the ground truth for estimating performance of the methods; and also for the policy risk mcasure Rpol 1 [E((|do(A)) 1), (xx)) 1) (),)) 1) E (Y]doA 0), π (x;,) 0){(5(z)) = 0)].
A brief overview the specifications of the models employed follows. All GPs employ RBF base kernel (also DKL's specifications in the last hidden layer). DKL and CounterDKL deep NN structure features three [10, 5, 2] hidden layers. The AutoEncoder deep structure employed for the "AutoEnc + GP and "AutoEnc + CounterGP' models similarly learns a 2-dimensional encoded lower-dimensional representation, where the encoder has two [10, 5] hidden layers before the 2-dim representation and the decoder has [5, 10] hidden layers before the reconstruction loss.
D Marginal Likelihood Maximization in Multioutput Deep Kernels
In the multitask deep kernel learning class of models, the parameter space è = (W, 6, B) is made of the deep neural network's weights W, the base kernel's hyperparameters á (variance, lengthscales, etc.) and the coregionalization matrix B entries. These parameters are learnt jointly by maximizing the log- marginall likelihood at the end of the GP layer. Using the chain rule, the derivatives are:
OC OL OK 0g(a, W)
aw oK, 0g(z, W)
ow
Oc OL 0Kg
0 0K, dΦ
OC OL aK
0B
OK
ÖB
where g( W) is the function mapping the inputs to the lower representation space parametrized by W, K is the base kernel and K(-) = ks(-) the coregionalizod kernel.
E Additional Simulated Experiments
Finally, we describe and present results on a few additional simulated examples that we conducted to assess CounterDKL performance compared to some other specifications seen in Section 5.2, on datasets with varying sample size, predictor space and action space dimensions. In particular, following Dudffk al. (2011) and Farajtabar al. (2018), we make use of some of the popular datasets for classification in the open-source UCT Machine Learning Repository (https://archive ics. uci..uci..oo www.m//indemnee.mm////// php), , by transforming the
19

---

Published in Transactions on Machine Learning Research (11/2022)
classification task in causal Off-Policy Evaluation task in the following way. Each dataset is equipped with pair of covariates X; and classification labels I4. We view the classification labels L; as our discrete actions L; = Ai, and consequently generate the action-specific outcome Yo, as function of the covariates as follows:
Ya. exp ( + t where (.,)
and Bn is a P-dimensional vector of action-specific coelficients, where entries are 0.. 0.0} sampled from Multinomial(O.5) 0.25, 0.15), with replacement. The datasets utilized are summarized in Table 3 in terms of sample size N, number of covariates P and number actions. We compare GP, CounterGP, DKL and CounterDKL models on an Off-Policy Evaluation task, where we evaluate the uniformly at randorn generated policy, via the absoulte regret or risk measure, defined as T.[ V(=e) - (TT.)) []. All models employ RBF base kernel, either directly the inputs replications of the experiments for each dataset are gathered in Table 4
the lower dimensional layer. Results averaged over B 20
20


Reviews:
[{"review_id": "tBFPVEj0Jq", "reviewer": "Reviewer_9sDy", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n - Combines three idea from the literature (CounterGP, DKL, and more general multi-task GP's) in a novel way to attack causal inference problems that are high-dimensional along multiple axes (many treatment levels, many outcomes, many covariates).\n - Experiment has a nice ablation design.\n\nWeaknesses:\n - In the presentation, too much time is spent reviewing background material. The key idea of using DKL is not introduced until page 7. Some of the background material is not particularly necessary or relevant (for example, there is no need to restate the backdoor adjustment formula as a theorem, and the re-statement of the minimax bound in Corollary 3.1 does not really speak to finite-sample performance---see requested changes).\n - Relatedly, the main technical contributions are difficult to parse out from the exposition. See requested changes for proposed focus.\n\n - The tradeoffs introduced by the modeling assumptions built into the framework are not well-described. There is a bias-variance tradeoff that comes from sharing information between surfaces (either between outcomes, or between counterfactuals). This is particularly relevant for sharing information between counterfactuals, because the correlation between counterfactuals is fundamentally unidentifiable from the data.\n - Relatedly, the experiments could be more thorough in exploring tradeoffs, particularly whether bias introduced by the GP prior really does fade with sample size.\n\n - Despite using a fully Bayesian framework and plotting credible bands in some of the demo figures, there is no evaluation of Bayesian uncertainty from the proposed method.\n - The experiments are not particularly compelling for the multi-outcome approach. It would be useful to at least specify a simulation where one can see the lift from the multi-outcome approach.\n - The authors also don't give an example of a real dataset that matches the multi-valued treatment multidimensional outcome setting that they built their method for.\n\n - There are several misstatements about causal inference with backdoor adjustment.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes a method for estimating counterfactual outcome surfaces using Gaussian processes that are parameterized using the DKL framework. The Gaussian processes are linked using a Kronecker structured covariance kernel that allows information sharing across counterfactuals that correspond to different treatments levels, and across potentially multi-dimensional outcomes. The authors compare the overall Multi-Output Counter DKL approach to ablated baselines in a simulation study, and compare the single outcome approach to previous GP-based approaches on the well-established LaLonde job training data.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "# Highlighting the main contributions\n\nThe main technical contributions are (1) introducing a new covariance structure that scales up the Kronecker product corresponding to the number of counterfactual outcome surfaces, and to adds a new Kronecker product corresponding to the outcome dimension, and (2) using the DKL structure to reduce the dimension of the GP on the covariate space, perhaps making the expanded covariance structure more tractable. I think focusing more of the paper on precisely describing the scaling challenges that this introduces and how DKL addresses them would be more appropriate. As of now, these technical contributions are difficult to pull out, where more than half of the paper is dedicated to restating known results.\n\nSlimming down the background material add adding more empirical exploration could be a useful way to reframe the paper.\n\nSome examples of unnecessary exposition:\n - Figure 1(b) and 1(c) are unnecessary. The working assumption is that $\\vX$ satisfies the backdoor criterion, so discussing other cases is not relevant to this paper.\n - Theorem 2.2 is not necessary (it is also wrong, see below). Simply quoting the backdoor adjustment formula would be sufficient.\n - Corollary 3.1 does not seem necessary. It is an asymptotic result, not finite sample, and the main point that it becomes hard to estimate an outcome surface when there are sparsely sampled input regions for a particular counterfactual does not require a formal statement. In fact, the target of estimation in this paper is not causal effects, but the intervention distributions themselves, so the primary content of this corollary (which relates the difficulty of estimating a contrast between surfaces to the difficulty of estimating each individual surface) is not relevant to the problem being solved. In terms of the challenges that result from many levels of A, simply stating that the sample size gets sharded across these actions, often very unevenly, seems to make this point effectively.\n\n# Exploring bias-variance tradeoffs\n \nThe authors note that the primary benefit of the multi-task GP framework is sharing information across outcome surfaces (either between counterfactuals or between outcome dimensions). However, this information sharing introduces bias. Right now the paper does not do a good job of noting the kind of risk that this bias can introduce, and the settings under which the bias would be beneficial. For example, in Figure 2, the red surface could have just as easily continued horizontally, in which case the multi-task GP estimate would be very far off.\n\nOne additional risk here is that the residual correlation between counterfactual outcomes---say, Y_a1, Y_a2---after conditioning on X is not identifiable, even with infinite data, but the proposed approach (and CounterGP) specify a correlation between them. If this specification introduces bias in the estimated intervened expectation surfaces, it would not fade away with larger sample size. This would also certainly affect the joint uncertainty bands associated with Y_a1 and Y_a2.\n\nIt would be useful to specify some simulations where the multi-task GP helps and some where it hurts (at least in finite samples), and see how the benefit fades with larger sample size.\n\n# More detailed empirical exploration\n\nThe empirical exploration in this paper could be much stronger. Evaluation of Bayesian uncertainty from the model, and evaluations on a wider range of simulations that showcase different aspects and tradeoffs of model performance would be a good start. Also, specifying a problem where the multi-outcome mutli-task GP clearly helps would be useful. Finally, demonstrating the method on a real dataset that has the multi-action-level multi-outcome structure would provide far stronger motivation. The authors claim that this is a common setting, so it should be doable to find such a dataset.\n\n# Misstatements in Section 2 on Causal Identification Conditions\n\nA DAG is not equivalent to an SCM. The DAG does not specify the F or p(u) part of the tuple. A better statement would be that an SCM implies a DAG.\n\nThe statement of Theorem 2.2 is incorrect. The backdoor adjustment formula says that p(\\vY | do(A=a), \\vX) = p(\\vY | A=a, \\vX). From this P(\\vY | do(A=a)) is recovered by E[p(\\vY | A=a, X)], where this expectation is taken over the marginal distribution of X, p(X). Note that P(\\vY | A=a) = E[p(\\VY | A=a, X) | A=a], which is very different, even if the backdoor criterion is satisfied. The last equality in the proof is wrong (it would only hold if p(x | A=a) = p(x)).\n\nTheorem 2.2 is also unnecessary. The statement of the do calculus puts no restrictions on the dimensionality of the variables in the graph, so this statement for multivariate Y is already implied by well-established results.\n\n“Provided that covariates X (or a subset of them) satisfy the backdoor criterion”. Note the famous M-bias or M-structure example in which a subset of X satisfy the backdoor criterion, but adjusting for all of X does not. I.e., satisfaction of the backdoor criterion of a subset of variables does not imply that the superset satisfies the backdoor criterion.\n\n# Other Comments\n\nThe quantity Y_i|do(A_i=a_i) is usually written as a counterfactual Y_{a} or a potential outcome Y(a). do notation is reserved for distributions (hence its use as a conditioning variable).\n\nThe LHS of the definition equation right before (4) does not really parse. Should we be conditioning on f_Y, or on m(.) and K(.,.)? If conditioning on f_Y, it seems like f_Y(x*_j) is just a constant. Meanwhile, eq (4) is written entirely in terms of the hyperparameters m and K.\n\nThe indexing of the coregionalization matrix is extremely difficult to parse. In the end, since Q is fixed to be 1, it seems like it would make sense to just present the intrinsic coregionalization matrix here, and note that it can be generalized to the LMC.\n\nThe LHS’s of the math environment on page 5 don’t seem to be consistent. In the second equation, should the arguments of k(.,.) be f_d(\\vx_i) or should they just be \\vx_i? It seems like the latter, given than the RHS is written in terms of \\vx_i.", "broader_impact_concerns": "No concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "6QiKzYk15Q", "reviewer": "Reviewer_6xmF", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "First of all, I remark that I liked the manuscript and enjoyed while reading it. The motivation, presentation, notation and development of the main ideas and methodology is very clear and thorough to me. Overall, the quality of the paper is high in its current state, with perhaps a few details or points for improvement. \n\n**Strengths**\n\n- The paper adequately identifies the benefits of introducing multi-output GPs in the counterfactual learning problem for exploiting the existing correlations.\n- Clarity is excellent in the manuscript, even for readers non-familiar with the specific problem of counterfactual learning.\n- The work is novel, as it extends *causal multi-task GPs* (Alaa & van der Schaar NIPS 2017) for accepting a stacked corregionalization model, and also introduces deep kernel learning for treating with high-dimensional covariates.\n- Experiments show empirical results of the performance and how the multi-output flavour helps when considering correlations between counterfactual functions.\n\n**Weaknesses**\n\n- Presentation of tasks based on *policies* could be better/longer explained for a better understanding, as it does not seem completelty clear to me.\n- Experimental results could be slightly extended to match one of the particular results included in (Alaa & van der Schaar 2017,2018). Since this is the main reference in the SOTA for the work.\n- Theoretical results in Section 3.2 are of interest for the manuscript, but perhaps *Corollary 3.1* could be better explained or at least explain somewhere how Eq. (5) is obtained.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper considers the problem of learning multiple counterfactual functions from observational data for causal inference. In particular, the manuscript focuses on the high-dimensional scenario with multiple actions and correlated outcomes. This becomes particularly convenient for multi-task Gaussian processes (GPs), that are scaled up using *deep kernels*. In the paper, the methodology for introducing multioutput deep kernels in the counterfactual learning problem is presented, followed by results showing the benefits of modelling the correlations among such functions. ", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "Additionally, I want to remark that authors are aware of the limitations that introducing DKL (Wilson et al.  2016) might have for tuning parameters, which are deeper discussed in Ober et al. (2021). Here, I would appreciate if authors could address/answer some of the questions and points of change that wrote down while reading the manuscript.\n\n**Questions**\n\n - **Q1** --  In the Introduction, it is written *\"exploration of policies in the real-world is usually very costly and potentially harmful\"*. Why harmful? Could you elaborate on a particular example?\n - **Q2** -- How restricted is the model to accept discrete actions -- i.e. binary. Could the authors indicate what would be the main changes for considering other continuous actions or similar?\n - **Q3** -- Why first equation in Section 3.3. is not like Eq. (2)? where did A_i go?\n - **Q4** -- Is CounterGP one of the models considered in Alaa & van der Schaar (2017,2018)?\n - **Q5** -- Can correlations vanish due the representation learning task (i) in the 2nd paragraph of the intro) of CounterDKL?\n\n**Changes**\n\n- **C1** -- I would love to see more details on Section 3.3., particularly on how both correlation matrices B_y and B_A emerge from the formulation. That would be fantastic for the paper.\n- **C2** -- Experiments are nice! but somehow baselines are a bit disconnected from the SOTA to me. Could authors put which model corresponds to what paper?\n- **C3** -- Experiment in Section 5.2. is important as it is the application one. A bit more of explanation on the particular properties of the data and the target task would help.\n- **C4** -- The paper claims a lot the DKL for high-dimensional scenarios, but later experiments are somehow on data with ~ 10 covariates. Would it be possible to show the performance on a larger dataset or task? At least to perceive the type of scalability of the counterfactual model.\n- **C5** -- I also think it would be important to compare the method with one of the experiments in Alaa & Schaar 2017, 2018 if it is possible. \n- **C6** -- Perhaps, explaining better or introducing the *policy-based* metrics in Section 5.1 would help.\n\n**Minors**\n\n- $\\epsilon$ corresponding to $\\mathcal{U}$ exogeneous variables is a bit confusing. I see that it is because of the latent functions $u$ in the MOGP. But perhaps there is a better choice...\n- *guarantees* just before Section 2.", "broader_impact_concerns": "I have no concerns on this", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "6jov57H8xY", "reviewer": "Reviewer_65W5", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n- Very good Introduction section.\n- The paper is written clearly and easy to read.\n- Empirical experiments are thorough and show the merits of the proposed method.\n\nWeaknesses\n- Contributions of the paper are stated under the Related Works section. The convention is that it’s stated at the end of the Introduction section.\n- Methods proposed in (Shalit et al, 2017), (Yao et al, 2018), etc. are not non-parametric. They are parametric methods, parameterized by neural networks.\n- In description of the overlap assumption, it is stated that $p(A_i=a | X_i=x) \\in (\\alpha, 1-\\alpha)$ where $\\alpha \\in (0, 0.5]$. Why \"$0.5]$\"? Shouldn’t this be \"$1)$\"? We only need non-zero probability on the entire support.\n- In Section 3.1, it is not clear what $Q$ and $R_q$ are, and how they are related. It is also not clear what $K_q$ is, and how it is related to $K$.\n- There is little elaboration on the training process: What are the learnable parameters? What is the objective function (if any)? Is the training end-to-end?\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper addresses the problem of performing causal inference in observational data, in settings where we have a large number of covariates, multiple discrete actions, and multiple correlated outcomes of interest. The authors propose CounterDKL (i.e., counterfactual multi-task Deep Kernel Learning method), which first learns a lower dimensional representation space of the covariates, and then places a multi-task Gaussian process (GP) on top of the learned representation. The simulated experiments demonstrate the benefits of the proposed method.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "- Discuss explicitly and clearly, what are the contributions made on top of (Alaa and van der Schaar, 2017, 2018), besides employing the Deep Kernel Learning component. \n- Address the mentioned weaknesses in the previous box.\n", "broader_impact_concerns": "N/A\n", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Two out of three reviewers noted the updated version sufficiently addressed their comments. 

The third review is more critical, especially on the issue of bias-variance tradeoffs: "...poses a number of modeling assumptions that introduce bias (in particular, bias that does not go away even with infinite samples) in order to make estimation more sample efficient. However, there is little to no discussion in the paper about how users of the method should reason about these design choices. Unlike standard ML settings, one cannot empirically calibrate bias-variance tradeoffs on causal estimation problems because ground truth individual causal effects are not observable. Thus, it is necessary to provide users with a picture of when modeling choices will succeed and, critically, when they will fail so that users can reason about the suitability of the method for their application. Neither the method development nor the experiments provide adequate resolution on this core issue."

My understanding is that this issue has been acknowledged by the authors, but perhaps some recommendations could be added based on empirical results to guide practitioners. 

One reviewer raised concerns about the new contribution list. I agree with this and would suggest having only the contributions in the list, not the background/structure of the paper.

All reviewers note the incremental modelling changes compared to prior work. However, as the paper is technically solid/correct and could be of interest to the TMLR community, I suggest accepting it with minor revision.

Decision: Accept (poster)