Controlled interventions provide the most direct source of information for learning causal effects. In particular, a dose-response curve can be learned by varying the treatment level and observing the corresponding outcomes. However, interventions can be expensive and time-consuming. Observational data, where the treatment is not controlled by a known mechanism, is sometimes available. Under some strong assumptions, observational data allows for the estimation of dose-response curves. Estimating such curves nonparametrically is hard: sample sizes for controlled interventions may be small, while in the observational case a large number of measured confounders may need to be marginalized. In this paper, we introduce a hierarchical Gaussian process prior that constructs a distribution over the doseresponse curve by learning from observational data, and reshapes the distribution with a nonparametric affine transform learned from controlled interventions. This function composition from different sources is shown to speed-up learning, which we demonstrate with a thorough sensitivity analysis and an application to modeling the effect of therapy on cognitive skills of premature infants. We introduce a new solution to the problem of learning how an outcome variable Y varies under different levels of a control variable X that is manipulated. This is done by coupling different Gaussian process priors that combine observational and interventional data. The method outperforms estimates given by using only observational or only interventional data in a variety of scenarios and provides an alternative way of interpreting related methods in the design of computer experiments. Many problems in causal inference [15] consist of having a treatment variable X and and outcome Y , and estimating how Y varies as we control X at different levels. If we have data from a randomized controlled trial, where X and Y are not confounded, many standard modeling approaches can be used to learn the relationship between X and Y . If X and Y are measured in an observational study, the corresponding data can be used to estimate the association between X and Y , but this may not be the same as the causal relationship of these two variables because of possible confounders. To distinguish between the observational regime (where X is not controlled) and the interventional regime (where X is controlled), we adopt the causal graphical framework of [17] and [20] . In Figure 1 we illustrate the different regimes using causal graphical models. We will use p(· | ·) to denote (conditional) density or probability mass functions. In Figure 1 (a) we have the observational, or "natural," regime where common causes Z generate both treatment variable X and outcome variable Y . While the conditional distribution p(Y = x | X = x) can be learned from this data, this quantity is not the same as p(Y = y | do(X = x)): the latter notation, due to Pearl [17] , denotes a regime where X is not random, but a quantity set by an intervention performed by an external agent. The relation between these regimes comes from fundamental invariance assumptions: when X is intervened upon, "all other things are equal," and this invariance is reflected by the fact that the model in Figure 1 Our goal is to learn the relationship where X ≡ {x 1 , x 2 , . . . , x T } is a pre-defined set of treatment levels. We call the vector f (X ) ≡ [f (x 1 ); . . . ; f (x T )] the response curve for the "doses" X . Although the term "dose" is typically associated with the medical domain, we adopt here the term dose-response learning in its more general setup: estimating the causal effect of a treatment on an outcome across different (quantitative) levels of treatment. We assume the causal structure information is known, complementing approaches for causal network learning [20, 10] by tackling the quantitative and statistical side of causal prediction. In Section 2, we provide the basic notation of our setup. Section 3 describes our model family. Section 4 provides a thorough set of experiments assessing our approach, including sensitivity to model misspecification. We provide final conclusions in Section 5. The target estimand p(Y = y | do(X = x)) can be derived from the structural assumptions of Figure 1 (b) by standard conditioning and marginalization operations: Notice the important difference between the above and p(Y = y | X = x), which can be derived from the assumptions in Figure 1 (a) by marginalizing over p(Z = z | X = x) instead. The observational and interventional distributions can be very different. The above formula is sometimes known as the back-door adjustment [17] and it does not require measuring all common causes of treatment and outcome. It suffices that we measure variables Z that block all "back-door paths" between X and Y , a role played by Z O in Figure 1 (c). A formal description of which variables Z will validate (2) is given by [21, 17, 20] . We will assume that the selection of which variables Z to adjust for has been decided prior to our analysis, although in our experiments in Section 4 we will assess the behavior of our method under model misspecification. Our task is to estimate (1) nonparametrically given observational and experimental data, assuming that Z satisfies the back-door criteria. One possibility for estimating (1) from observational data is consistent under some general assumptions on f (·) and g(·, ·). Estimating g(·, ·) nonparametrically seems daunting, since Z can in principle be high-dimensional. However, as shown by [6] , under some conditions 2 the problem of estimatingf (·) nonparametrically via (3) is no harder than a one-dimensional nonparametric regression problem. There is however one main catch: while observational data can be used to choose the level of regularization forĝ(·), this is not likely to be an optimal choice forf (·) itself. Nevertheless, even if suboptimal smoothing is done, the use of nonparametric methods for estimating causal effects by back-door adjustment has been successful. For instance, [8] uses Bayesian classification and regression trees for this task. Although of practical use, there are shortcomings to this idea even under the assumption that Z provides a correct back-door adjustment. In particular, Bayesian measures of uncertainty should be interpreted with care: a fully Bayesian equivalent to (3) would require integrating over a model for p(Z) instead of the empirical distribution for Z in D obs ; evaluating a dose x might require combining many g(x, z (i) ) where the corresponding training measurements x (i) are far from x, resulting on possibly unreliable extrapolations with poorly calibrated credible intervals. While there are well established approaches to deal with this "lack of overlap" problem in binary treatments or linear responses [19, 9] , it is less clear what to do in the continuous case with nonlinear responses. In this paper, we focus on a setup where it is possible to collect interventional data such that treatments are controlled, but where sample sizes might be limited due to financial and time costs. This is related to design of computer experiments, where (cheap, but biased) computer simulations are combined with field experiments [2, 7] . The key idea of combining two sources of data is very generic, the value of new methods being on the design of adequate prior families. For instance, if computer simulations are noisy, it is may not be clear how uncertainty at that level should be modeled. We leverage knowledge of adjustment techniques for causal inference, so that it provides a partially automated recipe to transform observational data into informed priors. We leverage knowledge of the practical shortcomings of nonparametric adjustment (3) so that, unlike the biased but low variance setup of computer experiments, we try to improve the (theoretically) unbiased but possibly oversmooth structure of such estimators by introducing a layer of pointtwise affine transformations. Heterogeneous effects and stratification. One might ask why marginalize Z in (2), as it might be of greater interest to understand effects at the finer subpopulation levels conditioned on Z. In fact, (2) should be seen as the most general case, where conditioning on a subset of covariates (for instance, gender) will provide the possibly different average causal effect for each given strata (different levels of gender) marginalized over the remaining covariates. Randomized fine-grained effects might be hard to estimate and require stronger smoothing and extrapolation assumptions, but in principle they could be integrated with the approaches discussed here. In practice, in causal inference we are generally interested in marginal effects for some subpopulations where many covariates might not be practically measurable at decision time, and for the scientific purposes of understanding total effects [6] at different levels of granularity with weaker assumptions. The main idea is to first learn from observational data a Gaussian process over dose-response curves, then compose it with a nonlinear transformation biased toward the identity function. The fundamental innovation is the construction a non-stationary covariance function from observational data. Given an observational dataset D obs of size N , we fit a Gaussian process to learn a regression model of outcome Y on (uncontrolled) treatment X and covariates Z. A Gaussian likelihood for Y given X and Z is adopted, with conditional mean g(x, z) and variance σ 2 g . A Matérn 3 /2 covariance function with automatic relevance determination priors is given to g(·, ·), followed by marginal maximum likelihood to estimate σ 2 g and the covariance hyperparameters [13, 18] . This provides a posterior distribution over functions g(·, ·) in the input space of X and Z. We then define f obs (X ), x ∈ X , as 3 where set {g(x, z (i) )} is unknown. Uncertainty about f obs (·) comes from the joint predictive distribution of {g(x, z (i) )} learned from D obs , itself a Gaussian distribution with a T N × 1 mean vector µ g and a T N × T N covariance matrix, T ≡ |X |. Since (4) is a linear function of {g(x, z (i) )}, this implies f obs (X ) is also a (nonstationary) Gaussian process with mean µ obs ( The motivation for (4) is that µ obs is an estimator of the type (3), inheriting its desirable properties and caveats. The cost of computing the covariance matrix K obs of f obs (X ) is O(T 2 N 2 ), potentially expensive. In many practical applications, however, the size of X is not particularly large as it is a set of intervention points to be decided according to practical real-world constraints. In our simulations in Section 4, we chose T = |X | = 20. Approximating such covariance matrix, if necessary, is a future research topic. Assume int chosen by some pre-defined design in X ). We assign a prior to f (·) according to the model where N (m, V) is the multivariate normal distribution with mean m and covariance matrix V, is the elementwise product, a(·) is a vector which we call the distortion function, and b(·) the translation function. The role of the "elementwise affine" transform a f obs +b is to bias f toward f obs with uncertainty that varies depending on our uncertainty about f obs . The multiplicative component a f obs also induces a heavy-tail prior on f . In the Appendix, we discuss briefly the alternative of using the deep Gaussian process of [5] in our observational-interventional setup. We parameterize K a as follows. Every entry k a (x, x ) of K a , (x, x ) ∈ X × X , assumes the shape of a squared exponential kernel modified according to the smoothness and scale information obtained from D obs . First, define k a (x, x ) as where (λ a , σ h ) are hyperparameters, δ(·) is the delta function, v x is a rescaling of K obs (x, x) 1/2 ,x is a rescaling of X to the [0, 1] interval,ŷ x is a rescaling of µ obs (x) to the [0, 1] interval. More precisely, Equation (6) is designed to borrow information from the (estimated) smoothness of f (X ), by decreasing the correlation of the distortion factors a(x) and a(x ) as a function of the Euclidean distance between the 2D points (x, µ obs (x)) and (x , µ obs (x )), properly scaled. Hyperparameter σ a controls how this distance is weighted. (6) also captures information about the amplitude of the distortion signal, making it proportional to the ratios of the diagonal entries of K obs (X ). Hyperparameter λ a controls how this amplitude is globally adjusted. Nugget 10 −5 brings stability to the sampling of a(X ) within Markov chain Monte Carlo (MCMC) inference. Hyper-hyperpriors on λ a and σ a are set as log(λ a ) ∼ N (0, 0.5), log(σ a ) ∼ N (0, 0.1). That is, λ a follows a log-Normal distribution with median 1, approximately 90% of the mass below 2.5, and a long tail to the right. The implied distribution for a(x) where s x = 1 will have most of its mass within a factor of 10 from its median. The prior on σ a follows a similar shape, but with a narrower allocation of mass. Covariance matrix K b is defined in the same way, with its own hyperparameters λ b and σ b . Finally, the usual Jeffrey's prior for error variances is given to σ Top row: scatterplot of observational data, with true dose-response function in solid green, adjusted µ obs in dashed red, and the unadjusted Gaussian process regression of Y on X in dashed-and-circle magenta (which is a very badly biased estimate in this example); scatterplot in the middle shows interventional data, 20 dose levels uniformly spread in the support of the observational data and 10 outputs per level − notice that the sign of the association is the opposite of the observational regime; matrix K obs is depicted at the end, where the non-stationarity of the process is evident. Middle row: priors constructed on f obs (X ) and a(X ) with respective means; plot at the end corresponds to the implied prior on a f obs + b. Bottom row: the respective posteriors obtained by Gibbs sampling. Figure 2 shows an example of inference obtained from synthetic data, generated according to the protocol of Section 4. In this example, the observational relationship between X and Y has the opposite association of the true causal one, but after adjusting for 15 of the 25 confounders that generated the data (10 confounders are randomly ignored to mimic imperfect prior knowledge), a reasonable initial estimate for f (X ) is obtained. The combination with interventional data results in a much better fit, but imperfections still exist at the strongest levels of treatment: the green curve drops at x > 2 stronger than the expected posterior mean. This is due to having both a prior derived from observational data that got the wrong direction of the dose-response curve at x > 1.5, and being unlucky at drawing several higher than expected values in the interventional regime for x = 3. The model then shows its strength on capturing much of the structure of the true dose-response curve even under misspecified adjustments, but the example provides a warning that only so much can be done given unlucky draws from a small interventional dataset. 5 In our experiments, we infer posterior distributions by Gibbs sampling, alternating the sampling of latent variables f (X ), a(X ), b(X ) and hyperparameters λ a , σ a , λ b , σ b , σ 2 int , using slice sampling [16] for the hyperparameters. The meaning of the individual posterior distribution over f obs (X ) might also be of interest. In principle, this quantity is potentially identifiable by considering a joint model for (D obs , D int ): in this case, f obs (X ) learns the observational adjustment g(x, z)p(z) dz. This suggests that the posterior distribution for f obs (X ) will change little according to model (5), which is indeed observed in practice and illustrated by Figure 2 . Learning the hyperparameters for K obs could be done jointly with the remaining hyperparameters, but the cost per iteration would be high due to the update of K obs . The MCMC procedure for (5) is relatively inexpensive assuming that |X | is small. Learning the hyperparameters of K obs separately is a type of "modularization" of Bayesian inference [11] . As we mentioned in Section 2, it is sometimes desirable to learn dose-response curves conditioned on a few covariates S ⊂ Z of interest. In particular, in this paper we will consider the case of straightforward stratification: given a set S of discrete covariates assuming instantiations s, we have functions f s (X ) to be learned. Different estimation techniques can be used to borrow statistical strength across levels of S, both for f s (X ) and f s obs (X ). However, in our implementation, where we assume |S| is very small (a realistic case for many experimental designs), we construct independent priors for the different f s obs (X ) with independent affine transformations. Finally, in the Appendix we also consider simple active learning schemes [12] , as suggested by the fact that prior information already provides different estimates of uncertainty across X (Figure 2 ), which is sometimes dramatically nonstationary. Assessing causal inference algorithms requires fitting and predicting data generated by expensive randomized trials. Since this is typically unavailable, we will use simulated data where the truth is known. We divide our experiments in two types: first, one where we generate random dose-response functions, which allows us to control the difficulty of the problem in different directions; second, one where we start from a real world dataset and generate "realistic" dose-response curves from which simulated data can be given as input to the method. We generate studies where the observational sample has N = 1000 data points and |Z| = 25 confounders. Interventional data is generated at three different levels of sample size, M = 40, 100 and 200 where the intervention space X is evenly distributed within the range shown by the observational data, with |X | = 20. Covariates Z are generated from a zero-mean, unit variance Gaussian with correlation of 0.5 for all pairs. Treatment X is generated by first sampling a function f i (z i ) for every covariate from a Gaussian process, summing over 1 ≤ i ≤ 25 and adding Gaussian noise. Outcome Y is generated by first sampling linear coefficients and one intercept to weight the contribution of confounders Z, and then passing the linear combination through a quadratic function. The dose-response function of X on Y is generated as a polynomial, which is added to the contribution of Z and a Gaussian error. In this way, it is easy to obtain the dose-response function analytically. Besides varying M , we vary the setup in three other aspects: first, the dose-response is either a quadratic or cubic polynomial; second, the contribution of X is scaled to have its minimum and maximum value spam either 50% or 80% of the range of all other causes of Y , including the Gaussian noise (a spam of 50% already generates functions of modest impact to the total variability of Y ); third, the actual data given to the algorithm contains only 15 of the 25 confounders. We either discard 10 confounders uniformly at random (the Random setup), or remove the "top 10 strongest" confounders, as measured by how little confounding remains after adjusting for that single covariate alone (the Adversarial setup). In the interest of space, we provide a fully detailed description of the experimental setup in the Appendix. Code is also provided to regenerate our data and re-run all of these experiments 1 . 1 Code available at http://www.homepages.ucl.ac.uk/~ucgtrbd/code/obsint. 6 Table 1 : For each experiment, we have either quadratic (Q) or cubic (C) ground truth, with a signal range of 50% or 80%, and an interventional sample size of M = 40, 100 and 200. E i denotes the difference between competitor i and our method regarding mean error, see text for a description of competitors. L i denotes the difference between our method and competitor i regarding log-likelihood (differences greater than 10 are ignored, see text). That is, positive values indicate our method is better according to the corresponding criterion. All results are averages over 50 independent simulations, italics indicate statistically significant differences by a two-tailed t-test at level α = 0. Evaluation is done in two ways. First, by the normalized absolute difference between an estimatef (x) and the true f (x), averaged over X . The normalization is done by dividing the difference by the gap between the maximum and minimum true values of f (X ) within each simulated problem 2 . The second measure is the log density of each true f (x), averaged over x ∈ X , according to the inferred posterior distribution approximated as a Gaussian distribution, with mean and variance estimated by MCMC. We compare our method against: I. a variation of it where a and b are fixed at 1 and 0, so the only randomness is in f obs ; II. instead of an affine transformation, we set f (X ) = f obs (X ) + r(X ), where r is given a generic squared exponential Gaussian process prior, which is fit by marginal maximum likelihood; III. Gaussian process regression with squared exponential kernel applied to the interventional data only and hyperparameters fitted by marginal likelihood. The idea is that competitors I and II provide sensitivity analysis of whether our more specialized prior is adding value. In particular, competitor II would be closer to the traditional priors used in computer-aided experimental design [2] (but for our specialized K obs ). Results are shown in Table 1 , according to the two assessment criteria, using E for average absolute error, and L for average log-likelihood. Our method demonstrated robustness to varying degrees of unmeasured confounding. Compared to Competitor I, the mean obtained without any further affine transformation already provides a competitive estimator of f (X ), but this suffers when unmeasured confounding is stronger (Adversarial setup). Moreover, uncertainty estimates given by Competitor I tend to be overconfident. Competitor II does not make use of our special covariance function for the correction, and tends to be particularly weak against our method in lower interventional sample sizes. In the same line, our advantage over Competitor III starts stronger at M = 40 and diminishes as expected when M increases. Competitor III is particularly bad at lower signal-tonoise ratio problems, where sometimes it is overly confident that f (X ) is zero everywhere (hence, we ignore large likelihood discrepancies in our evaluation). This suggests that in order to learn specialized curves for particular subpopulations, where M will invariably be small, an end-to-end model for observational and interventional data might be essential. 2 Data is also normalized to a zero mean, unit variance according to the empirical mean and variance of the observational data, in order to reduce variability across studies.  Figure 3: An illustration of a problem generated from a model fitted to real data. That is, we generated data from "interventions" simulated from a model that was fitted to an actual study on premature infant development [3] , where the dose is the number of days that an infant is assigned to follow a development program and the outcome is an IQ test at age 3. (a) Posterior distribution for the stratum of infants whose mothers had up to some high school education, but no college. The red curve is the posterior mean of our method, and the blue curve the result of Gaussian process fit with interventional data only. (b) Posterior distributions for the infants whose mothers had (some) college education. (c) The combined strata. We consider an adaptation of the study analyzed by [8] . Targeted at premature infants with low birth weight, the Infant Health and Development Program (IHDP) was a study of the efficacy of "educational and family support services and pediatric follow-up offered during the first 3 years of life" [3] . The study originally randomized infants into those that received treatment and those that did not. The outcome variable was an IQ test applied when infants reached 3 years. Within those which received treatment, there was a range of number of days of treatment. That dose level was not randomized, and again we do not have ground truth for the dose-response curve. For our assessment, we fit a dose-response curve using Gaussian processes with Gaussian likelihood function and the back-door adjustment (3) on available covariates. We then use the model to generate independent synthetic "interventional data." Measured covariates include birth weight, sex, whether the mother smoked during pregnancy, among other factors detailed by [8, 3] . The Appendix goes in detail about the preprocessing, including R/MATLAB scripts to generate the data. The observational sample contained 347 individuals (corresponding only to those which were eligible for treatment and had no missing outcome variable) and 21 covariates. This sample included 243 infants whose mother attended (some) high school but not college, and 104 with at least some college. We generated 100 synthetic interventional datasets stratified by mother's education, (some) high-school vs. (some) college. 19 treatment levels were pre-selected, amounting to 0 to 450 days with increments of 25 days. All variables were standardized to zero mean and unit standard deviation according to the observational distribution per stratum. Two representative simulated studies are shown in Figure 3 , depicting dose-response curves which have modest evidence of non-linearity, and differ in range per stratum 3 . On average, our method improved over the fitting of a Gaussian process with squared exponential covariance function that was given interventional data only. According to the average normalized absolute differences, the improvement was 0.06, 0.07 and 0.08 for the high school, college and combined data, respectively (where error was reduced in 82%, 89% and 91% of the runs, respectively), each in which 10 interventional samples were simulated per treatment level per stratum. We introduced a simple, principled way of combining observational and interventional measurements and assessed its accuracy and robustness. In particular, we emphasized robustness to model misspecification and we performed sensitivity analysis to assess the importance of each individual component of our prior, contrasted to off-the-shelf solutions that can be found in related domains [2] . 3 We do not claim that these curves represent the true dose-response curves: confounders are very likely to exist, as the dose level was not decided at the beginning of the trial and is likely to have been changed "on the fly" as the infant responded. It is plausible that our covariates cannot reliably account for this feedback effect. We are aware that many practical problems remain. For instance, we have not discussed at all the important issue of sample selection bias, where volunteers for an interventional study might not come from the same p(Z) distribution as in the observational study. Worse, neither the observational nor the interventional data might come from the population in which we want to enforce a policy learned from the combined data. While these essential issues were ignored, our method can in principle be combined with ways of assessing and correcting for sample selection bias [1] . Moreover, if unmeasured confounding is too strong, one cannot expect to do well. Methods for sensitivity analysis of confounding assumptions [14] can be integrated with our framework. A more thorough analysis of active learning using our approach, particularly in the light of possible model misspecification, is needed as our results in the Appendix only superficially covers this aspect.
The setting in REF differs significantly from ours as we do not assume access to any interventional data.
Observational-Interventional Priors for Dose-Response Learning
NIPS
Mathematics Computer Science