{'fraction_non_alphanumeric': 0.04304706132376207, 'fraction_numerical': 0.012633120026446905, 'mean_word_length': 4.6033342154009, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Counterfactual Learning with Multioutput Deep Kernels

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.


Introduction
In the recent years, there has been a surge of attention towards the use of machine learning methods for causal inference under observational data. The desire for highly personalized decision-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 at 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 effects on both treated and non-treated individuals. Thus, in the attempt to answer counterfactual questions about policy interventions, such as "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 in 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 observational data, with a focus on high-dimensional settings characterized by a 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 settings 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).
• We introduce the class of CounterDKL for causal inference under multiple actions/outcomes, discussing their benefits over causal multitask GPs in high-dimensional settings (in terms of covariates, actions and outcomes).
Studies with multiple outcomes are quite 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, in 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).

Related Work
Many significant contributions in the literature on non-linear regression techniques for counterfactual/causal learning have been made in the recent years. These include works on heterogeneous (or individual) treatment effects estimation (Shalit et al., 2017;Künzel et al., 2017;Alaa & van der Schaar, 2017;2018;Yao et al., 2018;Wager & Athey, 2018;Nie et al., 2020), and on the related field of off-policy evaluation and learning (Dudík et al., 2014;Farajtabar et al., 2018;Kallus, 2018;Athey & Wager, 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) 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-linear 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 or over-fitting of the surfaces of interest. Multitask GPs and DKL both allow for information 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).

Problem Framework
In order to introduce the problems of causal effects identification and estimation, we borrow the notation from do-calculus and Structural Causal Models (SCM) (Pearl, 2009), and start by defining the latter: Definition 2.1 (SCM). A Structural Causal Model is a 4-tuple E, V, F, p(ε) consisting of (subscript j indicates a random element in the set): 1. E: denoting a set of exogenous variables, defined as variables determined outside of the model.
2. V: denoting a set of endogenous variables, defined as variables determined inside the model.
3. F : a set of functions f j ∈ F mapping each element ε j ∈ E and every parent variables of V j ∈ V -which we denote by pa(V j ) -to the endogenous variables 4. p(ε j ): a probability distribution over ε j ∈ U.
A SCM implies a causal Directed Acyclic Graph (DAG), where nodes in the graph depict variables, i.e. ε j and V j , while edges denote the functional causal relationships f j ∈ F . In the following paragraphs, we examine the problem of identifying and estimating these causal relationships f j ∈ F , under observational data, relying on SCMs and causal DAGs. Figure 1: Causal DAG depicting the setting of interest, where arrows represent causal relationships. Set of observable covariates X ∈ X satisfies the backdoor criterion, in that conditioning on them allows for identification of the causal effects of action A ∈ A on the outcomes Y ∈ Y.

Suppose we have access to an observational dataset
set of discrete actions 2 (or manipulative variables), and Y i ∈ R M 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 Section 5.2. The overarching goal is to identify and estimate the (average) effects of intervening on the manipulative variable A i , by setting it equal to some value a, on the outcomes Y i . We denote the main quantity of interest, the joint interventional probability distribution of all the outcomes conditional on A i = 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 equations: where: f Y (·) is a vector-valued function if multiple outcomes are considered; pa(A i ) denotes parent variables (or causes) of A i ; ε i,j are error terms with a distribution p(ε i,j ). The equivalent causal DAG is represented in Figure 1, 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, or equivalently that the set of observed covariates X i ∈ X is causally sufficient, in the sense that conditioning on X i allows to identify the causal association between A i and Y i . Using Pearl's terminology, we equivalently say that X satisfies 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 known as positivity or common support assumption). Overlap requires that 0 < p(A i = a|X i = x) < 1, ∀a ∈ A -i.e. that the observed actions allocation given X i = x is never deterministic, so we could theoretically observe data points for which X i = 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 i in each action arm, so we can approximate extrapolation of model's prediction in those regions; thus, one must be careful as to which subpopulation, defined by a common support X over ⊆ X , to target to estimate causal effects. Under these two assumptions, the multivariate interventional distribution p(Y | do(A = a), X = x) can be recovered via backdoor adjustment (Pearl, 2009) (see Appendix A). Hence, provided that the covariates X i , direct common causes of A i and Y i , satisfy the backdoor criterion, we can estimate unbiased causal effects with the observed quantities in

Counterfactual Learning with Multitask GPs
For ease of exposition, let us consider the simple case depicted by the causal DAG in Figure 1a, with a single continuous outcome Y i ∈ R, with i ∈ {1, ..., N }. We tackle the problem of estimating p(Y |do(A = a)) via non-linear regression-adjustment (Johansson et al., 2016;Shalit et al., 2017;Künzel et al., 2017;Nie & Wager, 2020;Caron et al., 2022a). In particular, we assume, in line with most of the previous works, additive noise structure 3 , such that the outcome functional can be written as: There are different ways in which one can derive an estimator for p(Y |do(A = a)) and its moments, e.g. E(Y |do(A = a)), 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 i ∈ {0, 1}, as the quantity τ ( The idea behind causal multitask GPs is to view the D = |A| interventional quantities Y a , where D is the number of discrete action arms, as the output from a vector-valued function f Y (·) : X → R D (plus noise), modelled with a GP prior: with mean m(x i ) ∈ R D and covariance/kernel function K(x i , x j ) ∈ R D × R D , given two P -dimensional input points x i , x j ∈ X for units i and j. Given the likelihood function as a multivariate Gaussian on the diagonal and y i ∈ R D an output point, the posterior predictive distribution for a train set covariate realization x i ∈ X , train set outcome realization y i ∈ R and a test set covariate realization x * j ∈ X is obtained as, assuming zero prior mean m(·) = 0 for simplicity: and where φ denotes the model parameters and f * Y (x * j ) the function evaluated at a test point x * j . 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 a multitask GP is what induces task-relatedness when fitting the multi-valued surface f Y (·).

The multitask kernel
The simplest specification for the multitask kernel matrix is given by the separable kernels structure, which assumes single entries in K(·, ·) to be of the form k a,a (x i , with action a ∈ A = {1, ..., D}. Here, k(x i , x j ) represents a base kernel (e.g. linear, squared exponential, Matérn, etc.) while b a,a = k A (a, a ) is the generic entry of the D × D coregionalization matrix B, which contains the parameters governing task-relatedness over the actions A. In the trivial case where b a,a = 0 we have that tasks a and a are uncorrelated, i.e. actions a and a are unrelated in the way they affect the outcome Y .
Published in Transactions on Machine Learning Research (11/2022) A slightly more general framework, which we adopt in this work, is given by the sum of separable kernels structure (Álvarez et al., 2012). This assumes that the single entry of K(·, ·) reads k a,a (x i , i.e. the sum of Q different coregionalization matrices B q with associated base kernel k q (·, ·). In compact matrix notation this translates into K(X, X ) = Q q=1 B q ⊗ K q (X, X ) for two different input matrices X, X ∈ 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 Ya (·)} D a=1 generates from Q ≤ D common underlying independent latent GP functions {u q (·)} Q q=1 , parametrized by their base kernel k q (·, ·), (Álvarez et al., 2012).
In terms of the form of the coregionalization matrices B q , with q ∈ {1, ..., Q}, we will follow the linear model of coregionalization (LMC), which assumes that each B q is equal to B q = L q L q , with single entries b q; (a,a ) = Rq r=1 α r a,q α r a ,q . R q represents the number of GP samples obtained from the same latent GP function q, u q (·). Thus, adopting the LMC for causal learning is equivalent to assuming that correlation in the {f Ya (·)} D a=1 action-specific functions, modelled through the multitask kernel K(·, ·), arises from Q different samples of R q GP functions with the same kernel k q (·, ·), drawn from Q ≤ D different independent latent GP processes {u q (·)} Q q=1 . To express this more compactly, a causal multitask GP model under the 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 k a,a (x i , x j ) = BK(x i , x i ), with unique base kernelK(·, ·). 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 |X |, and direct likelihood maximization methods face issues in over-parametrized models, although some solutions, such as variational methods (Titsias, 2009;Hensman et al., 2013), might be adopted for better scalability.

Why multitask counterfactual learning?
We know that asymptotically the best approach to estimate the causal quantities f Y (x i ) 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 sample cases this is rarely the best strategy (Hahn et 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.e. 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 a and b, defined under the SCM specified by (2) as We assume that all the f a , ∀a ∈ A belong to the Hölder ball class of α a -smooth functions H(α a ), with α a − 1 bounded derivatives and bounded in sup-norm by a constant C > 0. Considering an L 2 -norm loss function on , the difficulty of CATE estimation with a nonparametric model ψ ∈ Ψ 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 covariate space is X = [0, 1] P , and f a depends on P a covariates s.t. P a ≤ P , ∀a ∈ A. Define n a,b < N as the subsample identified by action a and b. If both f a ∈ H(α a ) and Published in Transactions on Machine Learning Research (11/2022) , then CATE optimal minimax rate for a L 2 -norm loss function is: where x ∨ y = max{x, y} and is asymptotical equivalence. The first terms on the RHS of 6 relates to the problem of CATE function approximation, while the second term to the degree of sparsity of the CATE. Thus, 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 imbalanced), 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.e. regions in X where we mainly observe data points with specific action A i = a and very few others. In addition to this, as shown by Alaa & van der Schaar (2018), a multitask GP prior can achieve the optimal minimax rate of Corollary 3.1 in its posterior consistency.
Hence, splitting the sample into n a subgroups and fitting independent models can be very sample inefficient in these settings (6). Multitask GPs can aid extrapolation in such cases of strong sample selection bias, by learning the correlated functions {f Ya (·)} D a=1 jointly as f Y (·). The first row plots in Figure 2 provide a very simple one-covariate example of how multitask learning addresses the issue of extrapolation and prediction in poor overlap regions. Fitting the two surfaces f Y 1 (x i ) and f Y 0 (x i ) (dashed lines) through separate GP regressions results in a 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 f Ya (x i ) 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 2, the sample efficiency gains are less significant; 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; however, in settings with strong sample selection bias, or settings with multiple discrete actions or action spaces that grow with the sample size, the issue remains relevant.

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 i , in addition to the one over A i . The way we formulate this is through an additional coregionalization matrix, in the same fashion as we did for the actions case. This extended version featuring stacked coregionalization, has a GP prior of the following form: where B Y is the M × M coregionalization matrix over the outcomes, B A the D × D coregionalization matrix over the actions andK(·, ·) is the base kernel. The vector-valued function f Ya (·) in this case includes all the single-valued functionals {f a,m (·)} D,M a,m . The extra multitask learning layer over the outcomes Y i 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 it is also conceptually sound, as the quantity of interest is indeed the joint interventional distribution p(Y |do(A = a), X = x), which accounts for and explicitly models correlation between the outcomes, rather than the collection of marginal distributions {p(Y m |do(A = a))} M m=1 , which leaves correlation unspecified. As we will address in the later section, although the extra layer defined by B Y allows for higher sample efficiency, it also poses some issues due to parameter proliferation and stability of the optimization problem in high dimensions. Overlap is guaranteed 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 is generated as an increasing function of X i (selection bias). In the top row simulation, the two underlying counterfactual surfaces f Y d (x i ) (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.

Counterfactual Multitask Deep Kernel Learning
Gaussian Processes regressions are known to scale poorly with high dimensions. Their typical computational cost amounts to O(n 3 ) for training points and O(n 2 ) for test points. Similarly, coregionalized GPs suffer from over-parametrization and instability in the optimization procedure as the number of inputs P and the number of discrete actions D increase. Deep Kernel Learning (DKL) was firstly introduced by Wilson et al. (2016) with the aim of combining the scalability of Deep Learning methods and the nonparametric Bayesian uncertainty quantification of GPs in tackling prediction tasks in high-dimensional settings. Given a base kernel k(·, ·) (e.g. linear, squared exponential, etc.), a DKL structure learns a functional f Ya (·) by passing the P inputs X i ∈ X through a deep architecture (a fully-connected feedforward neural network in our case), which maps them to a lower dimensional representation space via non-linear activation functions. The base kernel k(·, ·) is then applied in this lower dimensional representation space, k(h (l) , h (l) ), where h (l) is a neural network's hidden layer, constituting a final Gaussian Process layer (or an infinite basis functions representation 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, namelyK(x, (Bohn et al., 2019). Thus, the DKL architecture is end-to-end, fully-connected and learnt jointly: the P inputs are passed on to hidden neural nets layers where the last hidden layer before the GP layer typically maps them to a lower dimensional representation space (with e.g. two hidden units). This is what generates benefits in terms of scalability compared to a 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 advantage of DKL is that the deep architecture preceeding the GP layer can itself learn arbitrarily complex function, so the choice of a specific GP kernel becomes less cumbersome. For example, Wilson et al. (2016) show that DKL is more robust in recovering step functions, due to weaker smoothness assumptions compared to standard GP kernels.
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 overly-parametrized (we refer to Ober et al. (2021) for a more detailed discussion of the issue). The kernel k(·, ·) in the last GP layer of a DKL architecture can easily incorporate the separable kernel structure for multitask learning, in the same fashion as the class of causal GPs presented earlier.  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.
to induce correlation across the action-specific functions {f Ya (·)} D a=1 , under the name of Counterfactual DKL (CounterDKL), where a similar Intrinsic Coregionalization Model (ICM) in the same spirit of (5) is placed on the last hidden layer of the neural network h (l) = g l (W l h (l−1) + b l ), where (W l , b l ) are the last layer's weights and bias, such that In this case, the Kronecker product of the coregionalization matrix occurs in the last hidden layer, and features lower dimensional representations instead of the potentially large number of raw inputs. Similarly, we can induce coregionalization over the M outcomes by adding another level of coregionalization, with the kernel reading K(h (l) , h (l) Figure  3 graphically depicts a counterfactual multitask Deep Learning architecture, with fully-connected hidden layers, a 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 i } l i=1 and biases {b i } l i=1 ; ii) the base kernelK hyperparameters Φ, e.g. in the case of squared exponential kernel Φ includes lengthscales and variances parameters Φ = { , σ 2 }; iii) the coregionalization matrix B entries. Hence, the parameter space is the collection Θ

Experiments
The fundamental problem of causal inference is that the interventional quantity p(Y |do(A = a), X i = x i ) 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 a data generating process with three different tasks, and on a real-world example combining experimental and observational data. For the first simulated experiment, we construct the DGP such that the backdoor criterion holds for X i ∈ X

Simulated Example
We consider a simulated setting with D = 4 possible actions A = {0, 1, 2, 3} and M = 2 correlated outcomes Y = (Y 1 , Y 2 ) ∈ R 2 . Actions and outcomes are generated according to a policy π b (x i ) = p(A i = a|X i = x i ) and an outcome function f Y (x i ), both dependent on the covariates X i ∈ 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 f Y d (·) distinctly for each outcome and for each action (GP); ii) counterfactual multitask GP regression (Alaa & van der Schaar, 2017;2018), with coregionalization over A i only, meaning that we fit two separate models for each outcome, but a unique multi-valued function model for A i (CounterGP); iii) counterfactual multioutput GP regression, a unique model with coregionalization both over A i and Y i (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, 50, 2] units, the DKL equivalent of ii) (CounterDKL); vi) counterfactual multioutput DKL, the DKL equivalent of iii) (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, 2000, 2500}. In the second version we fix sample size to N = 1500 and study the behaviour of the estimators with increasing number of covariates P ∈ {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 i = a on outcome Y i , given a certain realization of the covariates space, X i = x, i.e. the estimation of ICE: . This is carried out using a 80% training set, and evaluated via RMSE on a 20% left-out test set.
• OPE: The second is Off-Policy Evaluation, which is concerned with quantifying how good a given alternative policy π e is, compared to the action allocation policy that generated the data (behavior policy, π β ). This is done by estimating the policy value, defined as the cumulative reward V(π e ) = E X ,A,Y i π e (a i |x i ) Y i |do(A i = a i ) originating from π e . In our case we pick the alternative policy  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, represented by the γ parameter: higher values of γ corresponds to higher probability of being assigned to one of the two action A i = {3, 4}, thus generating more arms imbalance.
• OPL: 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: π * p ∈ arg max πp∈Π V(π p ). 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 π * (x i ) = a, i.e. the action that generates the best outcome for them.
Since we are dealing with M = 2 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 , we replicate the experiment B = 100 times to obtain Monte Carlo averages and 95% confidence intervals for the metrics. Results are depicted in Figure 4. 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-parametrization, as we had to omit it from the study of increasing predictors due to failed convergence for P > 10. The advantages over standard DKL regression instead are entirely attributable to sample efficiency gains from multitask coregionalization in CounterDKL and MODKL, both in the increasing N and increasing P studies. In the case of increasing P , we 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, 1 or 2 regularizers) in the deep architecture part.
In addition, we run a 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 estimates (then averaged over actions and outcomes). This is defined as wheref a (x i ) low andf a (x i ) upp are the lower and upper bands of the 95% credible interval output by the model onf a (x i ), while f a (x i ) is the true individual counterfactual outcome corresponding to action A i = a. Given fixed N = 2000 and P = 20 and a similar data generating process as before, we introduce the parameter γ, which governs the level of confounding, or the degree of groups imbalance in terms of action allocation. Particularly, for increasing values of γ, we assign higher probability of choosing one of the two Results are gathered in Figure 5, where MODKL display higher performance both in terms of error and uncertainty quantification.

Real-World Example: Job Training Programs and Unemployment
We demonstrate the efficiency of CounterDKL also on a second experiment taken from Shalit et al. (2017), 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 a randomized and an observational subgroups, where the aim is to estimate the effects of participation on a 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 i ∈ {0, 1} denotes participation to the job training program. The original outcome Y i is earnings after the program, which is censored continuous (Y i = 0 for unemployed units). However, following Shalit et al. (2017), we construct a binary indicator Y i ∈ {0, 1} 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/categorical 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 a degenerate Dirichlet distribution, which makes the GP classification task efficiently faster and more scalable. The 7 covariates X i ∈ X in the study are the following: age, 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 where T and C are the number of treated 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 π( . Notice that we cannot measure performance on ICE directly as this is always unobservable in real-world scenarios; also, we restrict analysis of average causal/treatment effects on the treated group since we 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. (2017). 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 1, in terms of 70%-30% train and test set Mean Absolute Error (MAE) on ATT, Policy Risk R pol and average runtime, accompanied by 10-fold cross-validated 95% error intervals. In this example multitasking is induced only over the binary treatment, as we deal with just a single outcome of interest. As the results depict, by operating jointly via a unique loss function, CounterDKL is significantly more Published in Transactions on Machine Learning Research (11/2022) efficient than naively applying dimensionality reduction and fitting a multitask GP on a lower dimensional space as two separate steps. It also displays gains over CounterGP, thanks to its deep component that guarantees better computational time (in terms of runtime) and scalability, and is 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).  The experiment is a semi-simulated setup. It makes use of real-world data from the IHDP study, a 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 a 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 i ∈ {0, 1} and the 25 available covariates X i ∈ X , but simulates the two continuous potential outcomes (Y 1 , Y 0 ) ∈ R 2 , as described in the non-linear "Response Surface B" setting in Hill (2011 Table 2 refers to 1000 replication of the experiment on 80%-20% train-test split as in Alaa & van der Schaar (2017).

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 it 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 a new class of counterfactual models based on 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/outcomes 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).

C The Job Training Data
The Job Training data (LaLonde, 1986) are a popular case study in the causal inference literature. They comprise a portion of data pertaining to a randomized experiment and a portion of observational data, with the randomized experiment 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 on the randomized units only) for the Average Treatment Effect on the Treated group (ATT) as ATT = T −1 Te i=1 y i − C −1 C i=1 y i , where T and C are the number of treated and control units in the experimental data, and treat this as the ground truth for estimating performance of the methods; and also for the policy risk measure R pol = 1 − E Y |do(A i = 1), π(x i ) = 1 p(π(x i ) = 1) + E Y |do(A i = 0), π(x i ) = 0 p(π(x i ) = 0) .
A brief overview on 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, φ, 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-marginal likelihood L at the end of the GP layer. Using the chain rule, the derivatives are: where g(x, W ) is the function mapping the inputs to the lower representation space parametrized by W , K φ is the base kernel and K(·) = B ⊗ K φ (·) is the coregionalized 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 Dudík et al. (2011) and Farajtabar et al. (2018), we make use of some of the popular datasets for classification in the open-source UCI Machine Learning Repository (https://archive.ics.uci.edu/ml/index.php), by transforming the Published in Transactions on Machine Learning Research (11/2022) classification task in a causal Off-Policy Evaluation task in the following way. Each dataset is equipped with a pair of covariates X i and classification labels L i . We view the classification labels L i as our discrete actions L i = A i , and consequently generate the action-specific outcome Y ai as function of the covariates as follows: Y ai = exp{X i β a } + ε i , where N (0, 0.5) and β a is a P -dimensional vector of action-specific coefficients, where entries are {0.4, 0.2, 0.0} sampled from a Multinomial(0.6, 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 of actions. We compare GP, CounterGP, DKL and CounterDKL models on an Off-Policy Evaluation task, where we evaluate the uniformly at random generated policy, via the absoulte regret or risk measure, defined as E | V(π e ) −V(π e ) | . All models employ a RBF base kernel, either directly on the inputs or on the lower dimensional layer. Results averaged over B = 20 replications of the experiments for each dataset are gathered in Table 4 . . .