Bayesian Nonparametric Causal Inference: Information Rates and Learning Algorithms

I Introduction

The problem of estimating the individualized causal effect of a particular intervention from observational data is central in many application domains and research fields, including public health and healthcare , computational advertising , and social sciences . With the increasing availability of data in all these domains, machine learning algorithms can be used to obtain estimates of the effect of an intervention, an action, or a treatment on individuals given their features and traits. For instance, using observational electronic health record datahttps://www.healthit.gov/sites/default/files/briefs/, machine learning-based recommender system can learn the individual-level causal effects of treatments currently deployed in clinical practice and help clinicians refine their current treatment policies . There is a growing interest in using machine learning methods to infer the individualized causal effects of medical treatments; this interest manifests in recent initiatives such as STRATOS , which focuses on guiding observational medical research, in addition to various recent works on causal effect inference by the machine learning community .

The problem of estimating individual-level causal effects is usually formulated within the classical potential outcomes framework, developed by Neyman and Rubin . In this framework, every subject (individual) in the observational dataset possesses two “potential outcomes”: the subject’s outcome under the application of the treatment, and the subject’s outcome when no treatment is applied. The treatment effect is the difference between the two potential outcomes, but since we only observe the “factual” outcome for a specific treatment assignment, and never observe the corresponding “counterfactual” outcome, we never observe any samples of the true treatment effect in an observational dataset. This is what makes the problem of causal inference fundamentally different from standard supervised learning (regression). Moreover, the policy by which treatments are assigned to subjects induces a selection bias in the observational data, creating a discrepancy in the feature distributions for the treated and control patient groups, which makes the problem even harder. Many of the classical works on causal inference have focused on the simpler problem of estimating average treatment effects, where unbiased estimators based on propensity score weighting were developed to alleviate the impact of selection bias on the causal estimands (see and the references therein).

While more recent works have developed machine learning algorithms for estimating individualized treatment effects from observational data in the past few years , the inference machinery built in most of these works seem to be rather ad-hoc. The causal inference problem entails a richer set of modeling choices and decisions compared to that of the standard supervised learning (regression) problem, which includes deciding what model to use, how to model the treatment assignment variables in the observational data, and how to handle selection bias, etc. In order to properly address all these modeling choices, one needs to understand the fundamental limits of performance in causal effect estimation problems, and how different modeling choices impact the achievable performance.

In this paper, we establish the fundamental limits on the amount of information that a learning algorithm can gather about the causal effect of an intervention given an observational data sample. We also provide guidelines for building proper causal inference models that “do not leave any information on the table” because of poor modeling choices. A summary of our results is provided in the following Section.

II Summary of the Results

We address the individualized causal effect estimation problem on the basis of the Neyman-Rubin potential outcomes model . We focus on Bayesian nonparametric learning algorithms, as they are immune to model mis-specification, and can learn highly heterogeneous response functions that one would expect to encounter in datasets with medical or social outcomes . In Section IV, we introduce the notion of information rate as a measure for the quality of Bayesian nonparametric learning of the individualized causal effects. The information rate is defined in terms of a measure of the Kullback-Leibler divergence between the true and posterior distributions for the causal effect. In Theorem 1, we establish the equivalence between Bayesian information rates and frequentist estimation rate. In the rest of the paper, we characterize: (1) the optimal information rates that can be achieved by any Bayesian nonparametric learning algorithm, and (2) the nature of the priors that would give rise to “informationally optimal” Bayesian nonparametric causal inference procedure.

In Section V, we establish the fundamental limit on the information rate that can be achieved by any Bayesian causal inference procedure using an information-theoretic lower bound based on Fano’s method. The optimal information rate is a property of the function classes to which the potential outcomes belong, and is independent of the inference algorithm. We show that the optimal information rate for causal inference is governed by the “rougher” of the two potential outcomes functions. We also show that the optimal information rates for causal inference are insensitive to selection bias (Theorem 2).

In Section VI, we characterize the Bayesian priors that achieve the optimal rate. We show that the most common modeling choice adopted in the literature, which is to augment the treatment assignment variable to the feature space, leads to priors that are suboptimal in terms of the achievable rate (Theorem 3). We show that informationally optimal priors are ones that place a probability distribution over a vector-valued function space, where the function space has its smoothness matching the rougher of the two potential outcomes functions. Since the true smoothness parameter of the potential outcomes functions is generally unknown a priori, we propose a prior adaptation procedure, called the information-based empirical Bayes procedure, which optimizes the Bayesian prior by maximizing an information-theoretic criterion on the recovered causal effects rather than maximizing the marginal likelihood of the observed (factual) data.

We conclude the paper by building an information-optimal Bayesian causal inference algorithm that is based on our analysis. The inference procedure embeds the potential outcomes in a vector-valued reproducing kernel Hilbert space (vvRKHS), and uses a multi-task Gaussian process prior (with a Matérn kernel) over that space to infer the individualized causal effects. We show that for such a prior, the proposed information-based empirical Bayes method exhibits an insightful factual bias and counterfactual variance decomposition. Experiments conducted on a standard dataset that is used for benchmarking causal inference models show that our model significantly outperforms the state-of-the-art.

III Related Work

We conduct our analysis within the potential outcomes framework developed by Neyman and Rubin . The earliest works on estimating causal effects have focused on the problem of obtaining unbiased estimates for the average treatment effects using observational samples. The most common well-known estimator for the average causal effect of a treatment is the propensity score weighting estimator, which simply removes the bias introduced by selection bias by giving weights to different samples that are inversely proportional to their propensity scores . More recently, the machine learning community has also developed estimators for the average treatment effects that borrows ideas from representation learning, i.e. see for instance the work in . In this paper, we focus on the individual, rather than the average causal effect estimation problem.

To the best of our knowledge, non of the previous works have attempted to characterize the limits of learning causal effects in either the frequentist or Bayesian setups. Instead, most previous works on causal effect inference have focused on model development, and various algorithms have been recently developed for estimating individualized treatment effects from observational data, mostly based on either tree-based methods , or deep learning methods . Most of the models that were previously developed for estimating causal effects relied on regression models that treat the treatment assignment variables (i.e. whether or not the intervention was applied to the subject) as an extended dimension in the feature space. Examples of such models include Bayesian additive regression trees (BART) , causal forests , balanced counterfactual regression , causal multivariate additive regression splines (MARS) , propensity-dropout networks , or random forests . In all these methods, augmenting the treatment assignment variable to the feature space introduces a mismatch between the training and testing distribution (i.e. covariate shift induced by the selection bias ). The different methods followed different approaches for handling the selection bias: causal forests use estimates of the propensity score for deriving a tree splitting rule that attempts to balance the treated and control populations, propensity-dropout networks use larger dropout regularization for training points with very high or very low propensity scores, whereas balanced counterfactual regression uses deep neural networks to learn a balanced representation (i.e. a feature transformation) that tries to alleviate the effect of the selection bias. Bayesian methods, like BART, do not address selection bias since the Bayesian posterior naturally incorporates uncertainty in regions of poor overlap in the feature space. As we show later in Sections VI and VIII, our analysis and experimental results indicated that, by augmenting the treatment assignment variable to the feature space, all these methods achieve a suboptimal information rate.

Our analysis is related to a long strand of literature that studied frequentist (minimax) estimation rates, or posterior contraction rates in standard regression problems . In Theorem 2, we show that the optimal information rate for causal inference has the same form as the optimal minimax estimation rate obtained by Stone in for standard nonparametric regression problems, when the true regression function is set to be the rougher of the two potential outcomes functions. Our analysis for the achievable information rates for Gaussian process priors uses the results by van Zanten and van der Vaart in .

IV Bayesian Nonparametric Causal Inference from Observational Data

In this section, we provide a general description for the Neyman-Rubin causal model considered in this paper (Subsection IV-A), and present the Bayesian nonparametric inference framework under study (Subsection IV-B).

Condition 1 (unconfoundedness): Treatment assignment decisions are independent of the outcomes given the subject’s features, i.e. $(Y^{(0)}_{i},Y^{(1)}_{i})\!\perp\!\!\!\perp\omega_{i}\,|\,X_{i}$.

Condition 2 (overlap): Every subject has a non-zero chance of receiving the treatment, and treatment assignment decisions are non-deterministic, i.e. $0<\gamma(x)<1$.

IV-B Bayesian Nonparametric Causal Inference

Throughout this paper, we consider the following signal-in-white-noise random design regression model for the potential outcomes:

A Bayesian procedure for estimating the ITE function entails specifying a prior distribution $\Pi$ over the response surfaces $f_{1}(x)$ and $f_{0}(x)$, which in turn induces a prior over $T(x)$. The nonparametric nature of inference follows from the fact that $\Pi$ is a prior over functions, and hence the estimation problem involves an infinite-dimensional parameter space. For a given prior $\Pi$, the Bayesian inference procedure views the observational dataset $\mathcal{D}_{n}$ as being sampled according to the following generative model:

Since we are interested in estimating an underlying true ITE function $T(x)$, we will analyze the Bayesian causal inference procedure within the frequentist setup, which assumes that the subjects’ outcomes $\{Y^{(\omega_{i})}_{i}\}^{n}_{i=1}$ are generated according to the model in (3) for a given true (and fixed) regression functions $f_{0}(x)$ and $f_{1}(x)$. That is, in the next Subsection, we will assess the quality of a Bayesian inference procedure by quantifying the amount of information the posterior distribution $d\Pi_{n}(T\,|\,\mathcal{D}_{n})=d\Pi_{n}(f_{1}-f_{0}\,|\,\mathcal{D}_{n})$ has about the true ITE function $T$. This type of analysis is sometimes referred to as the “Frequentist-Bayes” analysis .

IV-C Information Rates

How much information about the true causal effect function $T(.)$ is conveyed in the posterior $d\Pi_{n}(T\,|\,\mathcal{D}_{n})$? A natural measure of the “informational quality” of a posterior $d\Pi_{n}(T\,|\,\mathcal{D}_{n})$ is the information-theoretic criterion due to Barron , which quantifies the quality of a posterior via the Kullback-Leibler (KL) divergence between the posterior and true distributions. In that sense, the quality (or informativeness) of the posterior $d\Pi_{n}(T\,|\,\mathcal{D}_{n})$ at a feature point $x$ is given by the KL divergence between the posterior distribution at $x$, $d\Pi_{n}(T(x)\,|\,\mathcal{D}_{n})$, and the true distribution of $(Y^{(1)}-Y^{(0)})\,|\,X=x$. The overall quality of a posterior is thus quantified by marginalizing the pointwise KL divergence over the feature space $\mathcal{X}$. For a prior $\Pi$, true responses $f_{0}$ and $f_{1}$, propensity function $\gamma$, and observational datasets of size $n$, the expected KL risk is:

which by the convexity of the KL divergence in its second argument, and using Jensen’s inequality, is bounded above by

V Optimal Information Rates for Bayesian Causal Inference

In this Section, we establish a fundamental limit on the information rate that can be achieved by any sequence of posteriors $d\Pi_{n}(T\,|\,\mathcal{D}_{n})$ for a given causal inference problem. Let the achievable information rate for a given prior $\Pi$ and function classes $\mathcal{F}^{\alpha_{0}}$ and $\mathcal{F}^{\alpha_{1}}$, denoted by $I_{n}(\Pi;\mathcal{F}^{\alpha_{0}},\mathcal{F}^{\alpha_{1}},\gamma)$, be the rate obtained by taking the supremum of the information rate over functions in $\mathcal{F}^{\alpha_{0}}$ and $\mathcal{F}^{\alpha_{1}}$. This is a quantity that depends only on the prior but not on the specific realizations of $f_{0}$ and $f_{1}$. The optimal information rate is defined to be the maximum worst case achievable information rate for all functions in $\mathcal{F}^{\alpha_{0}}$ and $\mathcal{F}^{\alpha_{1}}$, and is denote by $I^{*}_{n}(\mathcal{F}^{\alpha_{0}},\mathcal{F}^{\alpha_{1}},\gamma)$. While the information rate $\mathcal{I}_{n}(\Pi;f_{0},f_{1},\gamma)$ characterizes a particular instance of a causal inference problem with $(f_{0},f_{1},\gamma)$ and a given Bayesian prior $\Pi$, the optimal information rate $I^{*}_{n}(\mathcal{F}^{\alpha_{0}},\mathcal{F}^{\alpha_{1}},\gamma)$ is an abstract (prior-independent) measure of the “information capacity” or the “hardness” of a class of causal inference problems (corresponding to response surfaces in $\mathcal{F}^{\alpha_{0}}$ and $\mathcal{F}^{\alpha_{1}}$). Intuitively, one expects that the limit on the achievable information rate will be higher for smooth (regular) response surfaces and for propensity functions that are close to 0.5 everywhere in $\mathcal{X}$. Theorem 2 provides a detailed characterization for the optimal information rates in general function spaces. Whether or not the Bayesian inference procedure achieves the optimal information rate will depend on the prior $\Pi$. In the next Section, we will investigate different design choices for the prior $\Pi$, and characterize the “capacity-achieving” priors that achieve the optimal information rate.

In Theorem 2, we will use the notion of metric entropy $H(\delta;\mathcal{F}^{\alpha})$ to characterize the “size” of general (nonparametric or parametric) function classes. The metric entropy $H(\delta;\mathcal{F}^{\alpha})$ of a function space $\mathcal{F}^{\alpha}$ is given by the logarithm of the covering number $N(\delta,\mathcal{F}^{\alpha},\rho)$ of that space with respect to a metric $\rho$, i.e. $H(\delta;\mathcal{F}^{\alpha})=\log(N(\delta,\mathcal{F}^{\alpha},\rho))$. A formal definition for covering numbers is provided below. Definition 2. (Covering number) A $\delta$-cover of a given function space $\mathcal{F}^{\alpha}$ with respect to a metric $\rho$ is a set of functions $\{f^{1},.\,.\,.,f^{N}\}$ such that for any function $f\in\mathcal{F}^{\alpha}$, there exists some $v\in\{1,.\,.\,.,N\}$ such that $\rho(f,f^{v})\leq\delta$. The $\delta$-covering number of $\mathcal{F}^{\alpha}$ is

The general characterization of the optimal information rates in Theorem 2 is cast into specific forms by specifying the regularity classes $\mathcal{F}^{\alpha_{0}}$ and $\mathcal{F}^{\alpha_{1}}$. Table I demonstrates the optimal information rates for standard function classes, including analytic, smooth, Hölder [30, Section 6.4], Sobolev , Besov [30, Section 6.3], and Lipschitz functions . A rough description for the optimal information rates of all nonparametric function spaces ($\alpha$-smooth, Hölder, Sobolev, Besov, and Lipschitz) can be given as follows. If $f_{0}$ is $\alpha_{0}$-regular (e.g. $\alpha_{0}$-differentiable) and $f_{1}$ is $\alpha_{1}$-regular, then the optimal information rate for causal inference is

where $\asymp$ denotes asymptotic equivalence, i.e. in Bachmann-Landau notation, $g(x)\asymp f(x)$ if $g(x)=\Theta(f(x))$. That is, the regularity parameter of the rougher response surface, i.e. $\alpha_{0}\wedge\alpha_{1}$, dominates the rate by which any inference procedure can acquire information about the causal effect. This is because, if one of the two response surfaces is much more complex (rough) than the other (as it is the case in the depiction in Figure 1), then the ITE function $T(x)$ would naturally lie in a function space that is at least as complex as the one that contains the rough surface. Moreover, the best achievable information rate depends only on the smoothness of the response surfaces and the dimensionality of the feature space, and is independent of the selection bias. Due to the nonparametric nature of the estimation problem, the optimal information rate for causal inference gets exponentially slower as we add more dimensions to the feature space .

Note that in Theorem 2, we assumed that for the surfaces $f_{0}$ and $f_{1}$, all of the $d$ dimensions of $\mathcal{X}$ are relevant to the two response surfaces. Now assume that surfaces $f_{0}$ and $f_{1}$ have relevant feature dimensions in the sets $\mathcal{P}_{0}$ and $\mathcal{P}_{1}$, respectively, where $|\mathcal{P}_{\omega}|=p_{\omega}\leq d,\,\omega\in\{0,1\}$ , then

where $\mathcal{F}^{\alpha_{\omega}}_{\mathcal{P}_{\omega}}$ denotes the space of functions in $\mathcal{F}^{\alpha_{\omega}}$ for which the relevant dimensions are in $\mathcal{P}_{\omega}$. In (7), the rate is dominated by the more complex response surface, where “complexity” here is manifesting as a combination of the number of relevant dimensions and the smoothness of the response over the those dimensions. One implication of (7) is that the information rate can be bottle-necked by the smoother of the response surfaces $f_{0}$ and $f_{1}$, if such a response has more relevant dimensions in the feature spaceA more general characterization of the information rate would consider the case when the responses have different smoothness levels on each of the $d$-dimensions. Unfortunately, obtaining such a characterization is technically daunting.. More precisely, if $\alpha_{0}<\alpha_{1}$, then the information rate can still be bottle-necked by the smoother surface $f_{1}$ as long as $p_{1}>\frac{\alpha_{1}}{\alpha_{0}}\,p_{0}$.

Since the optimal (Bayesian) information rate is a lower bound on the (frequentist) minimax estimation rate (Theorem 1), we can directly compare the limits of estimation in the causal inference setting (established in Theorem 2) with that of the standard nonparametric regression setting. It is well known that the optimal minimax rate for estimating an $\alpha$-regular function is $\Theta(n^{-2\alpha/(2\alpha+d)})$; a classical result due to Stone . The result of Theorem 2 (and the tabulated results in Table I) asserts that the causal effect estimation problem is as hard as the problem of estimating the “rougher” of the two surfaces $f_{0}$ and $f_{1}$ in a standard regression setup.

The fact that selection bias does not impair the optimal information rate for causal inference is consistent with previous results on minimax-optimal kernel density estimation under selection bias or length bias . In these settings, selection bias did not affect the optimal minimax rate for density estimation, but the kernel bandwidth optimization strategies that achieve the optimal rate needed to account for selection bias . In Section VI, we show that the same holds for causal inference: in order to achieve the optimal information rate, the strategy for selecting the prior $\Pi$ needs to account for selection bias. This means that even though the optimal information rates in the causal inference and standard regression settings are similar, the optimal estimation strategies in both setups are different.

VI Rate-adaptive Bayesian Causal Inference

In Section V, we have established the optimal rates by which any Bayesian inference procedure can gather information about the causal effect of a treatment from observational data. In this Section, we investigate different strategies for selecting the prior $\Pi$, and study their corresponding achievable information rates. (An optimal prior $\Pi^{*}$ is one that achieves the optimal information rate $I^{*}_{n}$.) A strategy for selecting $\Pi$ comprises the following three modeling choices:

How to incorporate the treatment assignment variable $\omega$ in the prior $\Pi$?

What function (regularity) class should the prior $\Pi$ place a probability distribution over?

What should be the smoothness (regularity) parameter of the selected function class?

The first modeling decision involves two possible choices. The first choice is to give no special role to the treatment assignment indicator $\omega$, and build a model that treats it in a manner similar to all other features by augmenting it to the feature space $\mathcal{X}$. This leads to models of the form

We refer to priors over models of the form above as Type-I priors. The second modeling choice is to let $\omega$ index two different models for the two response surfaces. This leads to models of the form ${\bf f}(x)=[f_{0}(x),f_{1}(x)]^{T},$ where $f_{0}\in\mathcal{F}^{\beta_{0}}$ and $f_{1}\in\mathcal{F}^{\beta_{1}}$ for some $\beta_{0},\,\beta_{1}>0$. We refer to priors over models of the form ${\bf f}(.)$ as Type-II priors.

Type I and II priors induce different estimators for $T(x)$. The posterior mean estimator for a Type-I prior is given by

The main difference between Type-I and Type-II priors is that the former restricts the smoothness of $f(x,\omega)$ on any feature dimension to be the same for $\omega=0$ and $\omega=1$. This also entails that the relevant dimensions for the two response surfaces ($\omega=0$ and $\omega=1$) need to be the same under a Type-I prior. (This is a direct consequence of the fact that Type-I priors give no special role to the variable $\omega$.) As a result, a priori knowledge (or even data-driven knowledge) on the differences between responses $f_{0}$ and $f_{1}$ (e.g. in terms of smoothness levels or relevant dimensions) cannot be incorporated in a Type-I prior. Type-II priors can incorporate such information as they provide separate models for $f_{0}$ and $f_{1}$. However, while Type-I priors give a posterior of $f_{0}$ and $f_{1}$ using a joint model that is fitted using all the observational data, Type-II priors use only the data for one population to compute posteriors of one response surface, which can be problematic if the two populations posses highly unbalanced relative sizes (e.g. treated populations are usually much smaller than control populations ).

Unlike their parametric counterparts, the nonparametric Type-I and II priors can (in general) learn the ITE function consistently, but how do their information rates compare? Subsection VI-A studies the achievable information rates for “oracle” Type-I and Type-II priors that are informed with the true smoothness parameters ($\alpha_{0}$ and $\alpha_{1}$) and relevant dimensions of the function classes $\mathcal{F}^{\alpha_{0}}$ and $\mathcal{F}^{\alpha_{1}}$. In Subsection VI-B, we study the (more realistic) setting when $\mathcal{F}^{\alpha_{0}}$ and $\mathcal{F}^{\alpha_{1}}$ are unknown, and investigate different strategies for adapting the prior $\Pi$ to the smoothness of the treated and control response surface in a data-driven fashion.

In this Subsection, we assume that the true smoothness and relevant dimensions for $f_{0}$ and $f_{1}$ are known a priori. In the following Theorem, we show that Type-II priors are generally a better modeling choice than Type-I priors. Theorem 3. (Sub-optimality of Type-I priors) Let $\Pi_{\beta}^{\circ}$ be the space of all Type-I priors that give probability one to draws from $H^{\beta}$, and let $\Pi_{\beta_{0,1}}^{\circ\circ}$ be the space of all Type-II priors that give probability one to draws from $(H^{\beta_{0}},H^{\beta_{1}})$. If $f_{0}\in H^{\alpha_{0}}_{\mathcal{P}_{0}}$, $f_{1}\in H^{\alpha_{1}}_{\mathcal{P}_{1}}$, and $\mathcal{P}_{0}\neq\mathcal{P}_{1}$, then

where $\gtrsim$ denotes asymptotic inequality. Proof. See Appendix B. ∎ Theorem 3 says that if $\mathcal{P}_{0}\neq\mathcal{P}_{1}$, then the information rate that any Type-I prior can achieve is always suboptimal, even if we know the relevant dimensions and the true smoothness of the response surfaces $f_{0}$ and $f_{1}$. The Theorem also says that an oracle Type-II prior can achieve the optimal information rate. When the the surfaces $f_{0}$ and $f_{1}$ have the same relevant dimensions and the same smoothness, the two priors achieve the same rate. More precisely, the best achievable information rate for a Type-I prior is given by

whereas for Type-II priors, the best achievable rate is

We note that most state-of-the-art causal inference algorithms, such as causal forests , Bayesian additive regression trees , and counterfactual regression , use Type-I regression structures for their estimates. The sub-optimality of Type-I priors, highlighted in Theorem 3, suggests that improved estimates can be achieved over state-of-the-art algorithms via a Type-II regression structure.

We now focus on the second and third modeling questions: on what function space should the prior be placed, and how should we set the regularity of the sample paths drawn from the prior? In the rest of this Section, we assume that the true response surfaces reside in Hölder spaces. One possible prior over Hölder balls is the Gaussian process $\mathcal{GP}(\mbox{Mat\'{e}rn}(\beta))$, with a Matérn covariance kernel and a smoothness parameter $\beta$. (Draws from such a prior are almost surely in a $\beta$-Hölder function space .) In the following Theorem, we characterize the information rates achieved by such a prior. Theorem 4. (The Matching Condition) Suppose that $f_{0}$ and $f_{1}$ are in Hölder spaces $H^{\alpha_{0}}$ and $H^{\alpha_{1}}$, respectively, and let

be a Type-II prior over $(H^{\beta_{0}},H^{\beta_{1}})$. If $(\beta_{0}\wedge\alpha_{0}\wedge\beta_{1}\wedge\alpha_{1})\geq d/2$, then we have that

where posterior consistency holds only if $\beta_{0}\leq\alpha_{0}$, and $\beta_{1}\leq\alpha_{1}$. Proof. See Appendix C. ∎ For a Type-I prior $\Pi(\beta)=\mathcal{GP}(\mbox{Mat\'{e}rn}(\beta))$, the upper bound on $I_{n}(\Pi(\beta);H^{\alpha_{0}},H^{\alpha_{1}})$ is $n^{\frac{-2\beta}{2\beta+d}}$, with consistency holding for $\beta\leq\alpha$. Using the results of the paper by Castillo in , the upper bound in Theorem 4 can be shown to be tight. Recall that the optimal information rate for causal inference in Hölder spaces is $I^{*}_{n}(H^{\alpha_{0}},H^{\alpha_{1}})=n^{\frac{-2(\alpha_{0}\wedge\alpha_{1})}{2(\alpha_{0}\wedge\alpha_{1})+d}}$ (Table I). Theorem 4 quantifies the information rates achieved by a Type-II prior with smoothness levels $\beta_{0}$ and $\beta_{1}$. The Theorem says that a prior can achieve the optimal information rate if and only if it captures the smoothness of the rougher of the two response surfaces. This gives rise to the following matching condition that a prior $\Pi(\beta_{0},\beta_{1})$ requires in order to provide an optimal rate:

That is, the regularity of the prior needs to match the rougher of the two surfaces, and the prior over the smoother surface needs to be at least as smooth as the rougher surface. Consistency holds only if the prior is at least as smooth as the true response, since otherwise the response surfaces would not be contained in the support of the prior. Note that Theorem 4 assumes that the true response surfaces exhibit a Hölder-type regularity, and that the prior $\Pi(\beta_{0},\beta_{1})$ is placed on a reproducing kernel Hilbert space with a particular kernel structure. While proving that the matching condition holds for general priors and function spaces is technically daunting, we believe that (given the results in Table I) the matching condition in Theorem 4 would hold for other notions of regularity (e.g. Sobolev, Lipschitz, etc), and for a wide range of practical priors. For instance, Theorem 4 holds for Gaussian processes with re-scaled squared exponential kernels .

To sum up this Subsection, we summarize the conclusions distilled from our analyses of the achievable rates for oracle priors. Priors of Type II are generally a better design choice compared to priors of Type I, especially when the two response surfaces exhibit different forms of heterogeneity. In order to achieve the optimal information rate, a typical condition is that the regularity of the prior needs to match that of the rougher of the two response surfaces. Since in practice we (generally) do not know the true smoothness of the response surfaces, we cannot build a prior that satisfies the matching condition. Practical causal inference thus requires adapting the prior to the smoothness of the true function in a data-driven fashion; we discuss this in the next Subsection.

VI-B Rate-adaptive Data-driven Priors

Note that, unlike in standard nonparametric regression, adapting the regularity of the prior for the causal inference inference task entails a mixed problem of testing and estimation, i.e. we need to test whether $\alpha_{0}$ is less than $\alpha_{1}$, and then estimate $\alpha_{0}$ (or $\alpha_{0}$). Hence, one would expect that the prior adaptation methods used in standard regression problems would not necessarily suffice in the causal inference setup. Prior adaptation can be implemented via hierarchical Bayes or empirical Bayes methods. Hierarchical Bayes methods specify a prior over $\beta=(\beta_{0},\beta_{1})$ (also known as the hyper-prior ), and then obtain a posterior over the regularity parameters in a fully Bayesian fashion. Empirical Bayes simply obtains a point estimate $\hat{\beta}_{n}$ of $\beta$, and then conducts inference via the prior specified by $\hat{\beta}_{n}$. We focus on empirical Bayes methods since the hierarchical methods are often impractically expensive in terms of memory and computational requirements. A prior $\Pi_{\hat{\beta}_{n}}$ induced by $\hat{\beta}_{n}$ (obtained via empirical Bayes) is called rate-adaptive if it achieves the optimal information rate, i.e. $I_{n}(\Pi_{\hat{\beta}_{n}})=I^{*}_{n}$.

In the rest of this Subsection, we show that marginal likelihood maximization, which is the dominant strategy for empirical Bayes adaptation in standard nonparametric regression , can fail to adapt to $\alpha_{0}\wedge\alpha_{1}$ in the general case when $\alpha_{0}\neq\alpha_{1}$. (This is crucial since in most practical problems of interest, the treated and control response surfaces have different levels of heterogeneity .) We then propose a novel information-based empirical Bayes strategy, and prove that it asymptotically satisfies the matching condition in Theorem 4. Finally, we conclude the Subsection by identifying candidate function spaces over which we can define the prior $\Pi$ such that we are able to both adapt to functions in Hölder spaces, and also conduct practical Bayesian inference in an algorithmically efficient manner.

The failure of likelihood-based empirical Bayes in the causal inference setup is not surprising as maximum likelihood adaptation is only optimal in the sense of minimizing the Kullback-Leibler loss for the individual potential outcomes. Optimal prior adaptation in our setup should be tailored to the causal inference task. Hence, we propose an information-based empirical Bayes scheme in which, instead of maximizing the marginal likelihood, we pick the smoothness level $\hat{\beta}_{n}$ that minimizes the posterior Bayesian KL divergence, i.e.

The information-based empirical Bayes estimator is simply a Bayesian estimator of $\beta$ with the loss function being the posterior KL risk in (4). Unlike the likelihood-based method, the objective in (8) is an direct measure for the quality of causal inference conducted with a prior $\Pi_{\beta}$. In the following Theorem, we show that $\hat{\beta}_{n}$ asymptotically satisfies the matching condition in Theorem 4. Theorem 5. (Asymptotic Matching) Suppose that $f_{0}$ and $f_{1}$ belong to the Hölder spaces $H^{\alpha_{0}}$ and $H^{\alpha_{1}}$, respectively, and let $\Pi(\beta)$ be a prior over Hölder space with order $\beta$. If $\hat{\beta}_{n}$ is obtained as in (8) using cross-validation, then under certain regularity conditions we have that $\hat{\beta}_{n}\overset{p}{\to}(\alpha_{0}\wedge\alpha_{1})$. Proof. See Appendix D. ∎ Theorem 5 says that the information-based empirical Bayes estimator is consistent. That is, the estimate $\hat{\beta}_{n}$ will eventually converge to $\alpha_{0}\wedge\alpha_{1}$ as $n\to\infty$. Note that this is a weaker result than adaptivity: consistency of $\hat{\beta}_{n}$ does not imply that the corresponding prior will necessarily achieve the optimal information rate. However, the consistency result in Theorem 5 is both strongly suggestive of adaptivity, and also indicative of the superiority of the information-based empirical Bayes method to the likelihood-based approach.

Note that, while the information-based empirical Bayes approach guarantees the asymptotic recovery of $\alpha_{0}\wedge\alpha_{1}$, it can still undersmooth the prior for the smoother response surface. This can problematic if we wish the posterior credible interval on $T(x)$ to be “honest”, i.e. possess frequentist coverage . A more flexible Type-II prior that assigns different smoothness parameters $\beta_{0}$ and $\beta_{1}$ to response surfaces $f_{0}$ and $f_{1}$ can potentially guarantee honest frequentist coverage in a manner similar to that provided by causal forests . As a consequence of Theorem 1, it turns out that the information-based empirical Bayes estimator in (8) is structurally similar to the risk-based empirical Bayes adaptation method proposed in . Hence, we conjecture that our proposed empirical Bayes procedure can guarantee frequentist coverage for the estimated causal effects under some conditions .

VI-B2 Concrete Priors for Bayesian Causal Inference

Assuming that $f_{0}$ and $f_{1}$ belong to Hölder spaces, what concrete priors should one use in order to achieve the optimal information rates? We have already shown (in Theorem 4) that the Gaussian process prior $\Pi(\beta)=\mathcal{GP}(\mbox{Mat\'{e}rn}(\beta))$, which places a probability distribution over a Hölder space with regularity $\beta$, can achieve the optimal information rate under the matching condition. Gaussian processes, in general, place a probability distribution on a reproducing kernel Hilbert space (RKHS) , the nature of which is determined by the kernel structure. It is worth mentioning that for kernels other than the Matérn kernel, the optimal rate might not be achievable. For instance, using the squared exponential kernel would lead to a suboptimal information rate of $(\log(n))^{-(\alpha_{0}\wedge\alpha_{1})/2+d/4}$, whereas spline kernels achieved a rate of $(n/\log(n))^{-2(\alpha_{0}\wedge\alpha_{1})/(2(\alpha_{0}\wedge\alpha_{1})+d)},$ which is optimal up to a logarithmic factor . In order for such kernels to achieve the optimal rates, their smoothness parameters (e.g. the length-scale parameter of the radial basis kernel) need to be re-scaled with the size of the observational data as explained in . Selection of the right kernel can be based either on prior knowledge of the response surfaces, or through a model selection procedure based on the information-theoretic criterion in (8).

Another possible option for Bayesian nonparametric priors place their probability mass on the space of piece-wise constant functions (trees) . The machine learning object operating on those spaces is the Bayesian additive regression trees (BART) algorithm, which was especially proven successful in causal inference problems, and was one of the winning algorithms in the 2016 Atlantic Causal Inference Conference Competitionhttp://jenniferhill7.wixsite.com/acic-2016. Since BART places a prior on a space of non-differentiable (piece-wise constant) functions, one would expect that the information rates achieved by BART would be inferior to those achieved by a GP. A carefully designed BART can only achieve a near-optimal rate of $(n/\log(n))^{-2(\alpha_{0}\wedge\alpha_{1})/(2\alpha_{0}\wedge\alpha_{1}+d)}$ . Our conclusion is that a Gaussian process is a better choice for causal modeling, not only because it can achieve better rates than BART, but also because its relatively tractable nature would allow for an easy implementation for the information-based empirical Bayes scheme in (8).

Finally, we note that if we know a priori which response surface is rougher, then prior adaptation can be achieved very easily by tuning the prior smoothness to the population that correspond to the rougher surface only. Such an adaptation can be done through the conventional marginal likelihood maximization method. It is worth mentioning though that while practitioners may know which surface is rougher a priori, it is less likely that in a high-dimensional space we would know ahead of time which variables are relevant to which surface. As we can see in the discussion after Theorem 3, the information rate is bottle-necked by complexity and not just smoothness. A smoother surface with more relevant dimensions can still bottle-neck the information rate. So practitioners should consider variable selection, and not just smoothness estimation, as a means to adapt the prior.

VII Practical Rate-adaptive Causal Inference with Multitask Gaussian Process Priors

The previous Section provided a detailed recipe for the informationally optimal Bayesian causal inference procedure. In particular, inference should be conducted through a Type-II Gaussian process prior on an RKHS space (Theorem 3 and Subsection VI-B). Moreover, the RKHS space should be defined through a Matérn covariance kernel with parameters $\beta_{0}$ and $\beta_{1}$ for response surfaces $f_{0}$ and $f_{1}$ (Subsection VI-B), and the parameters ${\bf\beta}=(\beta_{0},\beta_{1})$ should be optimized via the information-based empirical Bayes procedure in (8). In this Section, we construct a practical learning algorithm that follows this recipe.

We chose the Matérn covariance kernel as the underlying regularity of the vvRKHS since it can achieve the optimal information rate (see Appendix E). In order to avoid undersmoothing any of the two surfaces, we also chose to assign separate smoothness parameters $\beta_{0}$ and $\beta_{1}$ to $f_{0}$ and $f_{1}$, respectively. Standard intrinsic coregionalization models for vector-valued kernels impose the same covariance parameters for all outputs , which implies that the prior will have the same smoothness on both $f_{0}$ and $f_{1}$. Thus, we constructed a linear model of coregionalization (LMC) , which mixes two intrinsic coregionalization models as follows

where $k_{\omega}(x,x^{\prime})=\mbox{Mat\'{e}rn}(\beta_{\omega}),\,\omega\in\{0,1\},$ whereas ${\bf A}$ and ${\bf B}$ are given by

Now that we completely specified the multi-task GP prior for a given hyper-parameter set ${\bf\beta}$, the only remaining ingredient in the recipe is to implement the information-based empirical Bayes adaptation criterion in (8). The following Theorem gives an insightful decomposition of the information-based empirical Bayes objective for the multi-task GP model. (In the following Theorem, ${\bf Y^{(W)}}=[Y_{i}^{(\omega_{i})}]_{i}$ and ${\bf Y^{(1-W)}}=[Y_{i}^{(1-\omega_{i})}]_{i}$ are vectors comprising all factual and counterfactual outcomes associated with an observational dataset $\mathcal{D}_{n}$.) Theorem 6. (Factual bias and counterfactual variance decomposition) The minimizer ${\bf\beta}^{*}$ of the information-based empirical Bayes adaptation criterion in (8) is given by

where $\mbox{Var}_{\Pi_{\beta}}$ is the posterior variance and $\|.\|_{p}$ is the $p$-norm. Proof. See Appendix E. ∎ Theorem 6 states that, when the prior is specified as a multi-task GP, the information-based empirical Bayes criterion in (8) decomposes to factual bias and counterfactual variance termsThe objective function in Theorem 6 can be easily optimized via a leave-one-out cross-validation procedure. Refer to for a detailed explanation.. The factual bias term quantifies the empirical error in the observed factual outcome that results from selecting a particular smoothness level $\beta$. In that sense, the factual bias is a measure of the goodness-of-fit for the posterior mean resulting from a prior smoothness $\beta$. On the other hand, the counterfactual variance term quantifies the posterior uncertainty that would be induced in the unobserved counterfactual outcomes when selecting a smoothness level $\beta$. A small value for $\beta$ would lead to a rough posterior mean function, which corresponds to a good empirical fit for the data. On the contrary, a small value for $\beta$ would induce large uncertainty in the unobserved outcomes, which corresponds to large uncertainty in the counterfactual outcomes. The couterfactual variance thus acts as a regularizer for the factual bias that helps solving the joint testing-estimation problem of identifying the minimum of $\alpha_{0}$ and $\alpha_{1}$, and estimating the value of $\alpha_{0}\wedge\alpha_{1}$. That is, the regularizer attempts to protect the prior from falsely recognizing either $\alpha_{0}$ or $\alpha_{1}$ as being very low just because it over-fit the factual outcomes, and hence underestimating the true $\alpha_{0}\wedge\alpha_{1}$, thereby undersmoothing the prior and giving rise to a suboptimal information rate. The two terms work in opposite directions as shown in Figure 3(b): factual bias pushes for undersmoothed priors and counterfactual variance pushes for oversmoothed priors. Theorem 6 says that the resulting prior will lie on the optimal boundary in the large data limit.

Finally, we note that the factual bias and counterfactual variance trade-off automatically handles selection bias. That is, when there is a poor overlap between the treated and control populations, the posterior counterfactual variances would tend to be higher, and the information-based empirical Bayes method would tend to oversmooth the prior rather than fitting the factual data. Selection bias does not affect the optimal information rate, but it does affect the optimal strategy for achieving that rate as long as we decide to share parameters and data points between our models for the potential outcomes.

VIII Experiments

We sought to evaluate the finite-sample performance of the Bayesian causal inference procedure proposed in Section VII, and compare it with state-of-the-art causal inference models. Causal inference models are hard to evaluate , and obviously, it is impossible to validate a causal model using real-world data due to the absence of counterfactual outcomes. A common approach for evaluating causal models, which we follow in this paper, is to validate the model’s predictions/estimates in a semi-synthetic dataset for which artificial counterfactual outcomes are randomly generated via a predefined probabilistic model. To ensure a fair and objective comparison, we did not design the semi-synthetic dataset used in the experiments by ourselves, but rather used the (standard) semi-synthetic experimental setup designed by Hill in . In this setup, the features and treatment assignments are real but outcomes are simulated. The experimental setup was based on the IHDP dataset, a public dataset for data from a randomized clinical trial. We describe the dataset in more detail in the following Subsection.

The Infant Health and Development Program (IHDP) is an interventional program that is intended to enhance the cognitive and health status of low birth weight, premature infants through pediatric follow-ups and parent support groups . The semi-simulated dataset in is based on features for premature infants enrolled in a real randomized experiment that evaluated the impact of the IHDP on the subjects’ IQ scores at the age of three. Because the data was originally collected from a randomized trial, selection bias was introduced in the treatment assignment variable by removing a subset of the treated population. All outcomes (response surfaces) are simulated. The response surface data generation process was not designed to favor our method: we used the standard non-linear ”Response Surface B” setting in . The dataset comprises 747 subjects (608 control and 139 treated), and there are 25 features associated with each subject.

VIII-B Benchmarks

We compared our algorithm with various causal models and standard machine learning benchmarks which we list in what follows: $\clubsuit$ Tree-based methods (BART , causal forests (CF) , $\spadesuit$ Balancing counterfactual regression (balancing neural networks (BNN) , and counterfactual regression with Wasserstein distance metric (CFRW) ), $\bigstar$ Propensity-based and matching methods ($k$ nearest-neighbor ($k$NN), propensity score matching (PSM)), a $\diamondsuit$ nonparametric spline regression model (causal MARS ), and $\odot$ Doubly-robust methods (Targeted maximum likelihood (TML) ). We also compared the performance of our model with standard machine learning benchmarks, including linear regression (LR), random forests (RF), AdaBoost, XGBoost, and neural networks (NN). We evaluated two different variants of all the machine learning benchmarks: a $\square$ Type-I regression structure, in which we use the treatment assignment variable as an input feature to the machine leaning algorithm, and a $\otimes$ Type-II regression structure, in which we fit two separate models for treated and control populations. We compare all these benchmarks with our proposed model: a Type-II multi-task GP prior (MTGP) with a Matérn kernel optimized through information-based empirical Bayes. We also compare the proposed model with a Type-I multi-task GP model (with a Matérn kernel) optimized through likelihood-based empirical Bayes in order to verify the conclusions drawn from our analyses.

All machine learning benchmarks had their hyperparameters optimized via grid search using a held-out validation set. Hyper-parameter optimization was using the mean square error in the observed factual outcomes as the optimization objective. For BART, we used the default prior as in , and did not tune the model’s hyper-parameters. For BNN and CFRW, we used the neural network configurations reported in and . Causal MARS was implemented as described in . PSM was implemented as described in , and its performance was obtained by assuming that every patient’s estimated ITE is equal to the average treatment effect estimated by PSM. All benchmarks were implemented in Python, with the exception of BART, causal forests and TMLE, all of which were implemented in R. We used the R libraries bartMachine, grf, and tmle for the implementation of BART, causal forests and TMLE, respectively. Our method was implemented in Python using GPy, a library for Gaussian processes .

VIII-C Evaluation

VIII-D Results

As can be seen in Table II, the proposed Bayesian inference algorithm (Type-II MTGP) outperforms all other benchmarks in terms of the (in-sample and out-of-sample) PEHE. This result suggests that the proposed model was capable of adapting its prior to the data, and may have achieved the optimal (or a near-optimal) information rate. The PEHE results in Table II are the averages of 1000 experiments with 1000 different random realizations of the semi-synthetic outcome model. This means that our algorithm is consistently outperforming all other benchmarks as it is displaying a very tight confidence interval.

The benefit of the information-based empirical Bayes method manifests in the comparison with the Type-I MTGP prior optimized via likelihood-based empirical Bayes. The performance gain of the Type-II MTGP prior with respect to the Type-I MTGP prior results from the fact that the two response surfaces in the synthetic outcomes model have different levels of heterogeneity (the control response is non-linear whereas the treated response is linear. See the description of Response surface B in ). Our algorithm is also performing better than all other nonparametric tree-based algorithms. This is expected since, as we have discussed earlier in Subsection VI-B, an oracle BART prior can only achieve the optimal information rate up to a logarithmic factor. With the default prior, it is expected that BART would display a slow information rate as compared to our adapted, information-optimal Matérn kernel prior. Similar insights apply to the frequentist random forest algorithms, which approximates the true regression functions through non-differentiable, piecewise functions (trees), and hence is inevitably suboptimal in terms of the achievable minimax estimation rate.

Our model also outperforms all the standard machine learning benchmarks, whether the ones trained with a Type-I regression structure, or those trained with a Type-II structure. We believe that this is because our model outperforms the standard machine learning benchmarks since the information-based empirical Bayes method provides a natural protection against selection bias (via the counterfactual variance regularization). Selection bias introduces a mismatch between the training and testing datasets for all the machine learning benchmarks (i.e. a covariate shift ), and hence all machine learning methods exhibit high generalization errors.

IX Conclusions

In this paper, we studied the problem of estimating the causal effect of an intervention on individual subjects using observational data in the Bayesian nonparametric framework. We characterized the optimal Kullback-Leibler information rate that can be achieved by any learning procedure, and showed that it depends on the dimensionality of the feature space, and the smoothness of the “rougher” of the two potential outcomes. We characterized the priors that are capable of achieving the optimal information rates, and proposed a novel empirical Bayes procedure that is adapts the Bayesian prior to the causal effect function through an information-theoretic criterion. Finally, we used the conclusions drawn from our analysis and designed a practical Bayesian causal inference algorithm with a multi-task Gaussian process, and showed that it significantly outperforms the state-of-the-art causal inference models through experiments conducted on a standard semi-synthetic dataset.

Appendix A Proof of Theorem 2

Recall from (4) that the KL risk is given by

Using Pinsker’s inequality [31, Lemma 11.6.1], the KL divergence can be bounded below as follows:

where $\|.\|_{TV}$is the total variation distance between probability measures, which is given by the $L_{1}$ norm of the difference between $P(x)$ and $Q_{\mathcal{D}_{n}}(x)$ as follows:

Since the $L_{1}$ norm is bounded below by the $L_{2}$ norm, we can lower bound the KL divergence by combining (A.12) with Pinsker’s inequality as follows:

which leads to the following asymptotic inequality

Let $\delta_{\omega}$ be the solution to $H(\delta_{\omega};\,\mathcal{F}^{\alpha_{\omega}})\asymp n\,\delta^{2}_{\omega}$. We will prove that the optimal rate is $\Theta(\delta^{2}_{0}\vee\delta^{2}_{1})$ by first showing that $I^{*}_{n}(\mathcal{F}^{\alpha_{0}},\mathcal{F}^{\alpha_{1}})$ is lower bounded by, i.e. $I^{*}_{n}(\mathcal{F}^{\alpha_{0}},\mathcal{F}^{\alpha_{1}})=\Omega(\delta^{2}_{0}\vee\delta^{2}_{1}),$ and then show that $I^{*}_{n}(\mathcal{F}^{\alpha_{0}},\mathcal{F}^{\alpha_{1}})=O(\delta^{2}_{0}\vee\delta^{2}_{1})$. We start by observing that the causal inference problem can be described through the following Markov chain

The amount of information shared between the true function $T(.)$ and the estimate $\hat{T}(.)$ can be quantified by the mutual information $I(T;\hat{T})$. Given the Markov chain above, we can upper bound $I(T;\hat{T})$ as follows

where $(*)$ follows from the data processing inequality , and the supremum in $(\star)$ is taken over all possible priors. $I(T;\hat{T})$ is bounded below by the rate-distortion function

where ${\small(\bullet)}$ is an application of Markov’s inequality. By combining (A.20) with the result in (A.23), the lower bound in (A.18) can be further bounded below as follows

where $h(.)$ is the binary entropy. From (A.25), we have that

which is an incarnation of Fano’s inequality. By combining (A.24) with (A.26), we have the following inequality

From (A.23), the minimax risk $R^{*}_{\Pi}$ is bounded below by

Thus, the minimax risk can be bounded below as follows

Since $R^{*}_{\Pi}$ is strictly positive, then we have that

where $\delta$ is the solution to the transcendental equation

The metric entropy of a function space $\mathcal{F}^{\alpha_{\omega}}$ is given by $H(\delta,\mathcal{F}^{\alpha_{\omega}})=\log(N(\delta,\mathcal{F}^{\alpha_{\omega}})$, and hence (A.29) is written as

Since the metric entropy $H(\delta,\mathcal{F}^{\alpha_{\omega}})$ is a decreasing function of the smoothness parameter $\alpha_{\omega}$, then it follows that the solution $\delta^{*}$ of the transcendental equation in (A.30) is given by $\delta^{*}=\delta_{0}\vee\delta_{1}$, where $\delta_{\omega}$ is the solution to the equation

The equation in $(\ref{eqA14})$ has a solution for all $n$ when the function space $\mathcal{F}^{\alpha_{\omega}}$ has a polynomial or a logarithmic metric entropy , which is the case for all function spaces of interest (see Table I for evaluations of $\delta_{0}\vee\delta_{1}$ for various function spaces). It follows from (A.28) and (A.31) that

We now focus on upper bounding $R^{*}_{\Pi}$. From , we know that the minimax risk is upper bounded by the channel capacity in (A.16), which is further bounded above by the covering numbers as follows

For $\delta$ satisfying (A.31), we have that

and hence $R^{*}_{\Pi}\lesssim\delta^{2}_{0}\vee\delta^{2}_{1}$. It follows that

By combining (A.32) and (A.33), we have that $I^{*}_{n}=\Omega(\delta^{2}_{0}\wedge\delta^{2}_{1})$ and $I^{*}_{n}=O(\delta^{2}_{0}\vee\delta^{2}_{1})$, and hence it follows that

Appendix B Proof of Theorem 3

Note that when $f_{0}\in H^{\alpha_{0}}_{\mathcal{P}_{0}}$ and $f_{1}\in H^{\alpha_{1}}_{\mathcal{P}_{1}}$, the metric entropy of $H^{\alpha_{0}}_{\mathcal{P}_{0}}$ and $H^{\alpha_{1}}_{\mathcal{P}_{1}}$ are given by :

From Theorem 2, we know that the optimal information rate is given by $I^{*}_{n}(H^{\alpha_{0}}_{\mathcal{P}_{0}},H^{\alpha_{1}}_{\mathcal{P}_{1}})\asymp\delta_{0}^{2}\vee\delta_{1}^{2},$ where $\delta_{0}$ and $\delta_{1}$ are the solutions for $\delta^{\frac{-|\mathcal{P}_{\omega}|}{\alpha_{\omega}}}\asymp n\,\delta^{2}_{\omega},\,\omega\in\{0,1\}$. Thus, we have that

and hence the optimal information rate is given by

For a Type-II prior $\Pi\in\Pi_{\beta_{0,1}}^{\circ\circ}$ over the two Hölder spaces $H^{\alpha_{0}}_{\mathcal{P}_{0}}$ and $H^{\alpha_{1}}_{\mathcal{P}_{1}}$, with $\beta_{0}=\alpha_{0}$ and $\beta_{1}=\alpha_{1}$, the minimax estimation rates for nonparametric regression over $f_{0}$ and $f_{1}$ are

where the number of feature dimensions $|\mathcal{P}_{0}\cup\mathcal{P}_{1}|$ correspond to all the relevant dimensions for the regression function $f(x,\omega)$. The regression function on the discrete dimension $\omega$ can be estimated at the $\sqrt{n}$ parametric rate and hence it does not affect the minimax estimation rate given above. Since the rate $n^{\frac{-2(\alpha_{0}\wedge\alpha_{1})}{2\beta+|\mathcal{P}_{0}\cup\mathcal{P}_{1}|}}$ is strictly slower than the optimal rate of $n^{\frac{-2\alpha_{0}}{2\alpha_{0}+|\mathcal{P}_{0}|}}\vee n^{\frac{-2\alpha_{1}}{2\alpha_{1}+|\mathcal{P}_{1}|}}$ for all $\beta>0$, it follows that

Appendix C Proof of Theorem 4

Using Cauchy-Schwarz inequality, we obtain the following:

and similarly for $\|\sqrt{1-\gamma(x)}\cdot(\hat{f}_{0}(x)-f_{0}(x))\|^{2}_{2}$. The proof of the Lemma is concluded by observing that $\|\gamma(x)\|_{2}$ is $O(1)$ and $\|(\hat{f}_{1}(x)-f_{1}(x))^{2}\|_{2}\asymp\|(\hat{f}_{1}(x)-f_{1}(x))\|^{2}_{2}$. The same result can be arrived at via Minkowski inequality. ∎ Lemma C3. The support of the prior $\Pi(\beta)=\mathcal{GP}(\mbox{Mat\'{e}rn}(\beta))$ is the space of Hölder functions with order $\beta$. ∎ The proofs for Lemmas C1 and C3 are standard and can be found in and respectively.

From Lemma C2 and the equivalence in (C.37), we have that

where $\phi_{f_{\omega}}(\varepsilon)$ is the concentration function defined as :

The concentration function measures the amount of prior mass that $\Pi$ places around the true function $f_{\omega}$. The transcendental equation in (C.38) provides a valid contraction rate whenever consistency holds. Consistency of Bayesian inference holds whenever the true parameter (in this case the true function $f_{\omega}$) is in the support of the prior . From Lemmas C1 and C3, it follows that the necessary and sufficient conditions for consistency is that $\beta_{0}\leq\alpha_{0}$ and $\beta_{1}\leq\alpha_{1}$.

In [27, Lemma 4], the concentration function $\phi_{f_{\omega}}(\varepsilon)$ for a sufficiently smooth prior $\mathcal{GP}(\mbox{Mat\'{e}rn}(\beta_{\omega}))$, with $\beta_{\omega}>d/2$, and a sufficiently smooth true function $f_{\omega}\in H^{\alpha_{\omega}}$, with $\beta_{\omega}>d/2$, was obtained as follows:

Thus, combining (C.38) and (C.40), the posterior contraction rate for $\Pi(\beta_{\omega})$ around $f_{\omega}$ is the solution to:

which concludes the proof of the Theorem.

Appendix D Proof of Theorem 5

The empirical smoothness estimate $\hat{\beta}_{n}$ is obtained by minimizing the empirical objective:

The candidate in $\{\beta^{(1)},.\,.\,.,\beta^{(K_{n})}\}$ that minimizes the cross-validated risk $L^{k_{n}}(\mathcal{D}_{n})$ is

The consistency of the estimator $\beta^{(K^{*}_{n})}$ follows from the results of Dudoit and van der Laan on the asymptotic performance of model selection via cross-validation for general loss functions . Suppose that $\sup_{\mathcal{D}_{n},\beta}L(\mathcal{D}_{n},\beta)\leq\infty$, $\beta^{*}\in\{\beta^{(1)},.\,.\,.,\beta^{(K_{n})}\}$, and $\log(K_{n})/(\sqrt{nv}(L^{k_{n}}-L))\overset{p}{\to}0$ as $n\to\infty$. Then, from Theorem 2 in , we have that $L^{k_{n}}(\mathcal{D}_{n},\beta^{(k_{n})})-L(\mathcal{D}_{n},\beta)\overset{p}{\to}0$. Since $\beta^{*}=(\alpha_{0}\wedge\alpha_{1})$ is a unique minimizer of $L(\mathcal{D}_{n},\beta)$, then it follows from the argmin continuous mapping theorem for $M$-estimators that $\beta^{(K^{*}_{n})}\to(\alpha_{0}\wedge\alpha_{1})$ .

Appendix E Proof of Theorem 6

where the expectation in is taken with respect to ${\bf Y^{(1-W)}}|\mathcal{D}$. The Bayesian risk can be written as

The loss function $\hat{\mathcal{L}}$ conditional on a realization of the counterfactual outcomes is given by

References