Bayesian Nonparametric Causal Inference: Information Rates and Learning Algorithms
Ahmed M. Alaa Mihaela van der Schaar
stat.ME cs.LG stat.ML








mytempeqncnt
thmTheorem
lemLemma
prpProposition
crlryCorollary
defntnDefinition
illexampIllustrative example

cndtnCondition




proof[1][Proof]

  #1
definition[1][Definition]

  #1
example[1][Example]

  #1
remark[1][Remark]

  #1
















































  
  ./Figs/
  
  
  .pdf,.jpeg,.png

  
  
  
  
  
  
  
  






































arrows,calc


[2]()
  #1  #2


‖‖










































































































































op-tical net-works semi-conduc-tor





XXXXXXXX,  Vol. XX, No. X, XXXX 2018
Alaa et al.: Bayesian Nonparametric Causal Inference
Bayesian Nonparametric Causal Inference: Information Rates and Learning Algorithms
    Ahmed M. Alaa, Member, IEEE, and Mihaela van der Schaar, Fellow, IEEE

A. Alaa is with the Department of Electrical Engineering, University of California Los Angeles (UCLA), Los Angeles, CA, 90095, USA (e-mail: ahmedmalaa@ucla.edu). 
M. van der Schaar is with the Department of Engineering Science, University of Oxford, Parks Road, Oxford, OX1 3PJ, UK (e-mail: mihaela.vanderschaar@eng.ox.ac.uk). 
Section <ref> of this work was presented in part at the thirty-first annual conference on Neural Information Processing Systems (NIPS), 2017. This work was supported by the ONR Data science program.


    December 30, 2023
===========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================

 
We investigate the problem of estimating the causal effect of a treatment on individual subjects from observational data; this is a central problem in various application domains, including healthcare, social sciences, and online advertising. Within the Neyman-Rubin potential outcomes model, we use the Kullback-Leibler (KL) divergence between the estimated and true distributions as a measure of accuracy of the estimate, and we define the information rate of the Bayesian causal inference procedure as the (asymptotic equivalence class of the) expected value of the KL divergence between the estimated and true distributions as a function of the number of samples. Using Fano's method, we establish a fundamental limit on the information rate that can be achieved by any Bayesian estimator, and show that this fundamental limit is independent of the selection bias in the observational data. We characterize the Bayesian priors on the potential (factual and counterfactual) outcomes that achieve the optimal information rate.  As a consequence, we show that a particular class of priors that have been widely used in the causal inference literature cannot achieve the optimal information rate. On the other hand, a broader class of priors can achieve the optimal information rate.  We go on to propose a prior adaptation procedure (which we call the information-based empirical Bayes procedure) that 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. Building on our analysis, we construct an information-optimal Bayesian causal inference algorithm. This algorithm embeds the potential outcomes in a vector-valued reproducing kernel Hilbert space (vvRKHS), and uses a multi-task Gaussian process prior over that space to infer the individualized causal effects. We show that for such a prior, the proposed information-based empirical Bayes method adapts the smoothness of the multi-task Gaussian process to the true smoothness of the causal effect function by balancing a tradeoff between the factual bias and the counterfactual variance. We conduct experiments on a well-known real-world dataset and show that our model significantly outperforms the state-of-the-art causal inference models.     


Bayesian nonparametrics, causal effect inference, Gaussian processes, multitask learning, selection bias.  

 


§ 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 <cit.>, computational advertising <cit.>, and social sciences <cit.>. 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 data[https://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 <cit.>. 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 <cit.>, which focuses on guiding observational medical research, in addition to various recent works on causal effect inference by the machine learning community <cit.>.

The problem of estimating individual-level causal effects is usually formulated within the classical potential outcomes framework, developed by Neyman and Rubin <cit.>. 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 <cit.> 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 <cit.>, 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. 



§ SUMMARY OF THE RESULTS

We address the individualized causal effect estimation problem on the basis of the Neyman-Rubin potential outcomes model <cit.>. 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 <cit.>. In Section <ref>, 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 <ref>, 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 <ref>, 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.        



§ RELATED WORK


We conduct our analysis within the potential outcomes framework developed by Neyman and Rubin <cit.>. 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 <cit.>. 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 <cit.>. 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 <cit.>, or deep learning methods <cit.>. 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) <cit.>, causal forests <cit.>, balanced counterfactual regression <cit.>, causal multivariate additive regression splines (MARS) <cit.>, propensity-dropout networks <cit.>, or random forests <cit.>. 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 <cit.>). 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 <ref> and <ref>, 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 <cit.>. 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 <cit.> 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 <cit.>.    



§ 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 <ref>), and present the Bayesian nonparametric inference framework under study (Subsection <ref>). 





 §.§ The Neyman-Rubin Causal Model


Consider a population of subjects with each subject i possessing a d-dimensional feature X_i ∈𝒳. An intervention is applied to some subjects in the population: subject i's response to the intervention is a random variable denoted by Y^(1)_i, whereas the subject's natural response when no intervention is applied is denoted by Y^(0)_i. The two random variables, Y^(1)_i, Y^(0)_i ∈ℝ, are known as the potential outcomes. The causal effect of the intervention (treatment) on subject i is characterized through the difference between the two (random) potential outcomes (Y^(1)_i-Y^(0)_i) | X_i = x, and is generally assumed to be dependent on the subject's features X_i = x. Hence, we define the individualized treatment effect (ITE) for a subject i with a feature X_i = x as    

    T(x) = 𝔼[. Y^(1)_i-Y^(0)_i | X_i = x].

Our goal is to estimate the function T(x) from an observational dataset 𝒟_n, which comprises n independent samples of the random tuple {X_i, ω_i, Y^(ω_i)_i}, where ω_i ∈{0,1} is an intervention assignment indicator that indicates whether or not subject i has received the intervention (treatment) under consideration. The outcomes Y^(ω_i)_i and Y^(1-ω_i)_i are known in the literature as the factual and the counterfactual outcomes, respectively <cit.>, <cit.>. Intervention assignments generally depend on the subjects' features, i.e. ω_i ⊥⊥ X_i. This dependence is quantified via the conditional distribution ℙ(ω_i=1|X_i=x), also known as the propensity score of subject i <cit.>, <cit.>. In the rest of this paper, we denote the propensity score of a feature point x as γ(x).    

The observational dataset 𝒟_n = {X_i, ω_i, Y^(ω_i)_i}^n_i=1 is drawn from a joint density dℙ(X_i,ω_i,Y^(0)_i,Y^(1)_i), with a probability space (Ω,ℱ,ℙ) that supports the following standard conditions <cit.>:


  * Condition 1 (unconfoundedness): Treatment assignment decisions are independent of the outcomes given the subject's features, i.e. (Y^(0)_i,Y^(1)_i) ⊥⊥ω_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<γ(x)<1.
  


 §.§ Bayesian Nonparametric Causal Inference













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

    Y^(ω)_i = f_ω(X_i) + ϵ_i,ω, ω∈{0,1},

where ϵ_i,ω∼𝒩(0,σ^2_ω) is a Gaussian noise variable. It follows from (<ref>) that 𝔼[Y^(ω)_i |  X_i = x] = f_ω(x), and hence the ITE is given by T(x) = f_1(x)-f_0(x). The functions f_1(x) and f_0(x) correspond to the response surfaces over the subjects' feature space with and without the intervention; the difference between these two surfaces correspond to the individualized effect of the intervention. We assume that 𝒳 is a compact metric space (e.g. bounded, closed sets in ℝ^d), and that the true response surfaces f_ω: 𝒳→ℝ, ω∈{0,1} are totally bounded functions in a space of “smooth" or “regular" functions ℱ^α_ω, where α_ω is a smoothness (or regularity) parameter. This roughly means that f_ω is α_ω-differentiable; precise definitions for α_ω-regular function classes will be provided in subsequent Sections.

A Bayesian procedure for estimating the ITE function entails specifying a prior distribution Π 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 Π is a prior over functions, and hence the estimation problem involves an infinite-dimensional parameter space. For a given prior Π, the Bayesian inference procedure views the observational dataset 𝒟_n as being sampled according to the following generative model: 

    f_0, f_1    ∼Π,   X_i ∼ dℙ(X_i=x) 
    ω_i | X_i=x    ∼(γ(x)) 
    
    Y^(ω_i)_i | f_0, f_1,ω_i    ∼𝒩(f_ω_i(x), σ^2_ω_i),  i = 1,. . .,n.

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^(ω_i)_i}^n_i=1 are generated according to the model in (<ref>) 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Π_n(T | 𝒟_n) = dΠ_n(f_1-f_0 | 𝒟_n) has about the true ITE function T. This type of analysis is sometimes referred to as the “Frequentist-Bayes" analysis <cit.>.



 §.§ Information Rates

 
How much information about the true causal effect function T(.) is conveyed in the posterior dΠ_n(T | 𝒟_n)? A natural measure of the “informational quality" of a posterior dΠ_n(T | 𝒟_n) is the information-theoretic criterion due to Barron <cit.>, 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Π_n(T | 𝒟_n) at a feature point x is given by the KL divergence between the posterior distribution at x, dΠ_n(T(x) | 𝒟_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 𝒳. For a prior Π, true responses f_0 and f_1, propensity function γ, and observational datasets of size n, the expected KL risk is:     

    𝔻_n(Π; f_0, f_1,γ) = 𝔼_x[ 𝔼_𝒟_n[D_ KL(P(x)  Q_𝒟_n(x)).] ],

where D_KL(..) is the KL divergence[The KL divergence between probability measures P and Q is given by D_KL(PQ) = ∫log(dP/dQ) dP <cit.>. The existence of the Radon-Nikodym derivative dP/dQ_n in (<ref>) is guaranteed since P and Q are mutually absolutely continuous.], P(x) is the true distribution of T(x), i.e. P(x) = dℙ(Y^(1)-Y^(0) | X = x), and Q_𝒟_n(x) is the posterior distribution of T(x), given by:

    Q_𝒟_n(x)    = dΠ_n(Y^(1)-Y^(0) | X = x, 𝒟_n) 
        (⋆)= dΠ_n(T(x)+𝒩(0,σ^2_0+σ^2_1) | 𝒟_n), 
        (*)=∫𝒩(T(x),σ^2_0+σ^2_1)   dΠ_n(T(x) | 𝒟_n),

where steps (⋆) and (*) in (<ref>) follow from the sampling model in (<ref>). The expected KL risk 𝔻_n in (<ref>) marginalizes the pointwise KL divergence D_ KL(P(x)  Q_𝒟_n(x)) over the distribution of the observational dataset 𝒟_n (generated according to (<ref>)), and the feature distribution dℙ(X=x). Variants of the expected KL risk in (<ref>) have been widely used in the analysis of nonparametric regression models, usually in the form of (cumulative) Cesàro averages of the pointwise KL divergence at certain points in the feature space <cit.>. Assuming posterior consistency, the (asymptotic) dependence of 𝔻_n(Π; f_0, f_1, γ) on n reflects the rate by which the posterior dΠ_n(T | 𝒟_n) “sandwiches" the true ITE function T(x) everywhere in 𝒳. An efficient causal inference procedure would exhibit a rapidly decaying 𝔻_n: this motivates the definition of an “information rate".


   
Definition 1. (Information Rate) We say that the information rate of a Bayesian causal inference procedure is ℐ_n, for a sequence ℐ_n → 0, if 𝔻_n(Π; f_0, f_1, γ) is Θ(ℐ_n).   


   
Note that ℐ_n is the equivalence class of all sequences that have the same asymptotic rate of convergence. In the rest of this paper, we use the notation ℐ_n(Π; f_0, f_1,γ) to denote the information rate achieved by a prior Π in a causal inference problem instance described by the tuple (f_0, f_1,γ). The notion of an “information rate" for a Bayesian causal effect inference procedure is closely connected to the frequentist estimation rate (with respect to the L_2 loss) with T(.) as the estimand <cit.>. The following Theorem establishes such a connection.
    


Theorem 1. Let 𝔻_n(Π; f_0, f_1,γ) be the KL risk of a given Bayesian causal inference procedure, then we have that

    𝔻_n(Π; f_0, f_1,γ) ≤σ̅·𝔼_𝒟_n[ ‖𝔼_Π[ T | 𝒟_n ]-T‖^2_2 ],

for some σ̅ > 0, where ._2^2 is the L_2(ℙ)-norm with respect to the feature distribution, i.e. f_2^2 = ∫ f^2(x) dℙ(X=x).



Proof. Recall from (<ref>) that 𝔻_n is given by

    𝔻_n(Π; f_0, f_1,γ) = 𝔼_x[ 𝔼_𝒟_n[D_ KL(P(x)  Q_𝒟_n(x)).] ].

Based on (<ref>), D_ KL(P(x)  Q_𝒟_n(x).) can be written as

    D_ KL(P(x)  ∫𝒩(T(x),σ^2_0+σ^2_1)   dΠ_n(T(x) | 𝒟_n).),
 
which by the convexity of the KL divergence in its second argument, and using Jensen's inequality, is bounded above by

    D_ KL(P(x)  𝒩(𝔼_Π[ T(x) | 𝒟_n ],σ^2_0+σ^2_1).).

From the regression model in (<ref>), we know that P(x) = dℙ(Y^(1)-Y^(0) | X = x) ∼𝒩(T(x),σ^2_0+σ^2_1), and hence KL divergence above can be written as  

    D_ KL(𝒩(T(x),σ^2_0+σ^2_1)  𝒩(𝔼_Π[ T(x) | 𝒟_n ],σ^2_0+σ^2_1).),

which is given by 1/2(σ^2_0+σ^2_1)| 𝔼_Π[ T(x) | 𝒟_n ]-T(x) |^2 since D_ KL(𝒩(μ_0,σ^2) 𝒩(μ_1,σ^2)) = 1/σ^2|μ_1-μ_0|^2 <cit.>. Hence, the expected KL risk is bounded above as follows

    𝔻_n(Π; f_0, f_1,γ)    ≤𝔼_𝒟_n[ 𝔼_x[| 𝔼_Π[ T(x) | 𝒟_n ]-T(x) |^2/2(σ^2_0+σ^2_1)] ], 
       = 1/2(σ^2_0+σ^2_1) 𝔼_𝒟_n[ ‖ 𝔼_Π[ T | 𝒟_n ]-T ‖_2^2 ],

for all n ∈ℕ_+.   
 

 
Theorem 1 says that the information rate of causal inference lower bounds the rate of convergence of the L_2(ℙ) risk of the sequence of estimates T̂_n induced by the posterior mean ∫ T dΠ_n(T | 𝒟_n). The L_2(ℙ) risk ‖ 𝔼_Π[ T | 𝒟_n ]-T ‖_2^2 was dubbed the precision in estimating heterogeneous effects (PEHE) by Hill in <cit.>, and is the most commonly used metric for evaluating causal inference models <cit.>. Theorem 1 tells us that the PEHE is Ω(ℐ_n), and hence the Bayesian information rate presents a limit on the achievable performance of frequentist estimation. In that sense, the asymptotic behavior of ℐ_n(Π; f_0, f_1,γ) is revealing of both the informational quality of the Bayesian posterior, as well as the convergence rates of frequentist loss functions.



§ 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Π_n(T | 𝒟_n) for a given causal inference problem. Let the achievable information rate for a given prior Π and function classes ℱ^α_0 and ℱ^α_1, denoted by I_n(Π; ℱ^α_0, ℱ^α_1, γ), be the rate obtained by taking the supremum of the information rate over functions in ℱ^α_0 and ℱ^α_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 ℱ^α_0 and ℱ^α_1, and is denote by I^*_n(ℱ^α_0, ℱ^α_1, γ). While the information rate ℐ_n(Π; f_0, f_1, γ) characterizes a particular instance of a causal inference problem with (f_0, f_1, γ) and a given Bayesian prior Π, the optimal information rate I^*_n(ℱ^α_0, ℱ^α_1, γ) 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 ℱ^α_0 and ℱ^α_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 𝒳. 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 Π. In the next Section, we will investigate different design choices for the prior Π, and characterize the “capacity-achieving" priors that achieve the optimal information rate.   

In Theorem 2, we will use the notion of metric entropy H(δ; ℱ^α) to characterize the “size" of general (nonparametric or parametric) function classes. The metric entropy H(δ; ℱ^α) of a function space ℱ^α is given by the logarithm of the covering number N(δ,ℱ^α,ρ) of that space with respect to a metric ρ, i.e. H(δ; ℱ^α) = log(N(δ,ℱ^α,ρ)). A formal definition for covering numbers is provided below.


   
Definition 2. (Covering number) A δ-cover of a given function space ℱ^α with respect to a metric ρ is a set of functions {f^1,. . ., f^N} such that for any function f ∈ℱ^α, there exists some v ∈{1,. . .,N} such that ρ(f, f^v) ≤δ. The δ-covering number of ℱ^α is <cit.> 

    N(δ,ℱ^α,ρ) inf{N ∈ℕ: }.

That is, the covering number of a function class ℱ^α is the number of balls (in a given metric ρ) of a fixed radius δ > 0 required to cover it. Throughout this paper, the metric entropy will always be evaluated with respect to the L_2(ℙ) norm. In the light of the definition above, the metric entropy can be thought of as a measure of the complexity of a function class; smoother function classes would generally display a smaller metric entropy. All function classes considered in this paper have finite metric entropy. Figure <ref> shows a pictorial depiction for two exemplary function classes ℱ^α_0 and ℱ^α_1 for the treated and control responses, respectively. In this depiction, α_0 is smaller than α_1, hence the δ-cover of ℱ^α_0 contains more balls than the δ-cover of ℱ^α_1, and it follows that ℱ^α_0 has a larger metric entropy. This manifests in the control response surface f_0 being less smooth than the treated response surface f_1. This is usually the case for real-world data on responses to medical treatments, where the untreated population typically display more heterogeneity than the treated population <cit.>.



We now present the main result of this Section. In the following Theorem, we provide a general characterization for the optimal information rates of Bayesian causal inference when the treated and control surfaces are known to belong to function classes ℱ^α_1 and ℱ^α_0.



Theorem 2. (Optimal Information Rates) Suppose that 𝒳 is a compact subset of ℝ^d, and that Conditions 1-2 hold. Then the optimal information rate is Θ(δ^2_0 ∨δ^2_1), where δ_ω is the solution for H(δ_ω; ℱ^α_ω) ≍ n  δ^2_ω, ω∈{0,1}. 
  

 
Proof. See Appendix A.     



Theorem 2 characterizes I^*_n(ℱ^α_0, ℱ^α_1,γ) in terms of the metric entropies H(δ; ℱ^α_0) and H(δ; ℱ^α_1) for general function classes ℱ^α_0 and ℱ^α_1. We used the local Fano method to derive an information-theoretic lower bound on the information rate that can be achieved by any estimator <cit.>. The characterization in Theorem 2 implies that selection bias has no effect on the achievable information rate. (Thus, in the rest of the paper we drop the dependency on γ when referring to I^*_n.) That is, as long as the overlap condition holds, selection bias does not hinder the information rate that can be achieved by a Bayesian causal inference procedure, and we can hope to find a good prior Π that achieves the optimal rate of posterior contraction around the true ITE function T(x) irrespective of the amount of bias in the data. Theorem 2 also says that the achievable information rate is bottle-necked by the more “complex" of the two response surfaces f_0 and f_1. Hence, we cannot hope to learn the causal effect at a fast rate if either of the treated or the control response surfaces are rough, even when the other surface is smooth.  

The general characterization of the optimal information rates in Theorem 2 is cast into specific forms by specifying the regularity classes ℱ^α_0 and ℱ^α_1. Table <ref> demonstrates the optimal information rates for standard function classes, including analytic, smooth, Hölder <cit.>, Sobolev <cit.>, Besov <cit.>, and Lipschitz functions <cit.>. A rough description for the optimal information rates of all nonparametric function spaces (α-smooth, Hölder, Sobolev, Besov, and Lipschitz) can be given as follows. If f_0 is α_0-regular (e.g. α_0-differentiable) and f_1 is α_1-regular, then the optimal information rate for causal inference is 

    I^*_n(ℱ^α_0, ℱ^α_1) ≍ n^-2(α_0 ∧α_1)/2(α_0 ∧α_1)+d,
  
where ≍ denotes asymptotic equivalence, i.e. in Bachmann-Landau notation, g(x) ≍ f(x) if g(x) = Θ(f(x)). That is, the regularity parameter of the rougher response surface, i.e. α_0 ∧α_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 <ref>), 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 <cit.>. 

Note that in Theorem 2, we assumed that for the surfaces f_0 and f_1, all of the d dimensions of 𝒳 are relevant to the two response surfaces. Now assume that surfaces f_0 and f_1 have relevant feature dimensions in the sets 𝒫_0 and 𝒫_1, respectively, where |𝒫_ω|=p_ω≤ d, ω∈{0,1} <cit.>, then

    I^*_n(ℱ^α_0_𝒫_0, ℱ^α_1_𝒫_1) ≍ n^-2α_0/2α_0+p_0∨ n^-2α_1/2α_1+p_1,
  
where ℱ^α_ω_𝒫_ω denotes the space of functions in ℱ^α_ω for which the relevant dimensions are in 𝒫_ω. In (<ref>), 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 (<ref>) 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 space[A 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 α_0 < α_1, then the information rate can still be bottle-necked by the smoother surface f_1 as long as p_1 > α_1/α_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 α-regular function is Θ(n^-2α/(2α+d)); a classical result due to Stone <cit.>. The result of Theorem 2 (and the tabulated results in Table <ref>) 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 <cit.>. 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 <cit.>. In Section <ref>, we show that the same holds for causal inference: in order to achieve the optimal information rate, the strategy for selecting the prior Π 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. 














§ RATE-ADAPTIVE BAYESIAN CAUSAL INFERENCE


In Section <ref>, 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 Π, and study their corresponding achievable information rates. (An optimal prior Π^* is one that achieves the optimal information rate I^*_n.) A strategy for selecting Π comprises the following three modeling choices: 


  * How to incorporate the treatment assignment variable ω in the prior Π?  

  * What function (regularity) class should the prior Π 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 ω, and build a model that treats it in a manner similar to all other features by augmenting it to the feature space 𝒳. This leads to models of the form 

    f(x,ω): 𝒳×{0,1}→ℝ.

We refer to priors over models of the form above as Type-I priors. The second modeling choice is to let ω index two different models for the two response surfaces. This leads to models of the form f(x) = [f_0(x),f_1(x)]^T, where f_0 ∈ℱ^β_0 and f_1 ∈ℱ^β_1 for some β_0, β_1 > 0. We refer to priors over models of the form 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

    T̂_n(x) = 𝔼_Π[f(x,1) | 𝒟_n]-𝔼_Π[f(x,0) | 𝒟_n],

whereas for a Type-II prior, the posterior mean ITE estimator is given by T̂_n(x) = 𝔼_Π[ f^T(x)  v | 𝒟_n], where v = [-1, 1]^T. Figure <ref> is a pictorial depiction for the posterior mean ITE estimates obtained via Type-I and Type-II priors. 


The main difference between Type-I and Type-II priors is that the former restricts the smoothness of f(x,ω) on any feature dimension to be the same for ω = 0 and ω = 1. This also entails that the relevant dimensions for the two response surfaces (ω = 0 and ω = 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 ω.) 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 <cit.>).    

In order to better illustrate the difference between Type-I and Type-II priors, we look at their simpler parametric counterparts. A Type-I linear regression model defines f(x,ω) as a linear function Y = β^T  x + γ̃ · ω + ε, where β∈ℝ^d, γ̃∈ℝ, and ε is a Gaussian noise variable. (Here the Type-I prior is a prior on the model coefficients β and γ̃) As we can see, this model restricts the two responses f_0 and f_1 to have the exact same interactions with the features through the coefficients in β. If we know a priori that f_0 and f_1 have different “slopes" or different relevant dimensions, we cannot incorporate this knowledge into the model. What would such a model learn? Assuming consistency, the estimated ITE function would be T̂_n(x) →γ̃ everywhere in 𝒳. Thus, the restricted nature of a Type-I parametric model led to a constant (non-individualized) estimate of T(.). On the contrary, a Type-II model of the form Y^(ω) = β_ω^T  x + ε, ω∈{0,1}, would allow for learning a linear estimate T̂_n(x) of the true function T(x), with potentially different relevant dimensions for both surfaces. However, Type-II model will only use data with ω=w to fit the model for Y^(w),  w ∈{0,1}.    

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 <ref> studies the achievable information rates for “oracle" Type-I and Type-II priors that are informed with the true smoothness parameters (α_0 and α_1) and relevant dimensions of the function classes ℱ^α_0 and ℱ^α_1. In Subsection <ref>, we study the (more realistic) setting when ℱ^α_0 and ℱ^α_1 are unknown, and investigate different strategies for adapting the prior Π to the smoothness of the treated and control response surface in a data-driven fashion. 



 §.§ Oracle Priors
 

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 Π_β^∘ be the space of all Type-I priors that give probability one to draws from H^β, and let Π_β_0,1^∘∘ be the space of all Type-II priors that give probability one to draws from (H^β_0,H^β_1). If f_0 ∈ H^α_0_𝒫_0, f_1 ∈ H^α_1_𝒫_1, and 𝒫_0 ≠𝒫_1, then

    inf_β_0,β_1inf_Π∈Π_β_0,1^∘∘ I_n(Π;H^α_0_𝒫_0,H^α_1_𝒫_1) ≍ I^*_n(H^α_0_𝒫_0,H^α_1_𝒫_1),
    inf_βinf_Π∈Π_β^∘ I_n(Π;H^α_0_𝒫_0,H^α_1_𝒫_1) ≳ I^*_n(H^α_0_𝒫_0,H^α_1_𝒫_1),

where ≳ denotes asymptotic inequality.


 
Proof. See Appendix B.     


  
Theorem 3 says that if 𝒫_0 ≠𝒫_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    

    inf_βinf_Π∈Π_β^∘ I_n(Π;H^α_0_𝒫_0,H^α_1_𝒫_1) = Θ(n^-2(α_0 ∧α_1)/2(α_0 ∧α_1)+|𝒫_0 ∪𝒫_1|),
 
whereas for Type-II priors, the best achievable rate is 

    inf_β_0,β_1inf_Π∈Π_β_0,1^∘∘ I_n(Π;H^α_0_𝒫_0,H^α_1_𝒫_1) = Θ(n^-2α_0/2α_0+|𝒫_0|∨ n^-2α_1/2α_1+|𝒫_1|).
 
We note that most state-of-the-art causal inference algorithms, such as causal forests <cit.>, Bayesian additive regression trees <cit.>, and counterfactual regression <cit.>, 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 𝒢𝒫((β)), with a Matérn covariance kernel and a smoothness parameter β. (Draws from such a prior are almost surely in a β-Hölder function space <cit.>.) 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^α_0 and H^α_1, respectively, and let

    ,

be a Type-II prior over (H^β_0,H^β_1). If (β_0 ∧α_0 ∧β_1 ∧α_1) ≥ d/2, then we have that

    I_n(Π(β_0,β_1);H^α_0,H^α_1) ≲ n^-2β_0/2β_0 + d∨ n^-2β_1/2β_1 + d,
    
where posterior consistency holds only if β_0 ≤α_0, and β_1 ≤α_1.



Proof. See Appendix C.     


      
For a Type-I prior Π(β) = 𝒢𝒫((β)), the upper bound on I_n(Π(β);H^α_0,H^α_1) is n^-2β/2β + d, with consistency holding for β≤α. Using the results of the paper by Castillo in <cit.>, 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^α_0,H^α_1) =  n^-2(α_0 ∧α_1)/2(α_0 ∧α_1) + d (Table <ref>). Theorem 4 quantifies the information rates achieved by a Type-II prior with smoothness levels β_0 and β_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 Π(β_0,β_1) requires in order to provide an optimal rate:

    [1.1]β_ω = α_0 ∧α_1, α_ω≤β_1-ω≤α_1-ω, ω = min_w ∈{0,1}α_w.

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 Π(β_0,β_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 <cit.>.

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.   
   


 §.§ 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 α_0 is less than α_1, and then estimate α_0 (or α_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 β = (β_0,β_1) (also known as the hyper-prior <cit.>), and then obtain a posterior over the regularity parameters in a fully Bayesian fashion. Empirical Bayes simply obtains a point estimate β̂_n of β, and then conducts inference via the prior specified by β̂_n. We focus on empirical Bayes methods since the hierarchical methods are often impractically expensive in terms of memory and computational requirements. A prior Π_β̂_n induced by β̂_n (obtained via empirical Bayes) is called rate-adaptive if it achieves the optimal information rate, i.e. I_n(Π_β̂_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 <cit.>, can fail to adapt to α_0 ∧α_1 in the general case when α_0 ≠α_1. (This is crucial since in most practical problems of interest, the treated and control response surfaces have different levels of heterogeneity <cit.>.) 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 Π 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.
 



  §.§.§ Information-based Empirical Bayes


To see why the marginal likelihood-based empirical Bayes method may fail in adapting priors for causal inference, consider the following example. Suppose that f_0 ∈ H^α_0 and f_1 ∈ H^α_1, where α_0 < α_1. Let Π(β̂_n) be a Type-I data-driven prior, where β̂_n is an empirical Bayes estimate of α_0 ∧α_1. For the likelihood-based empirical Bayes, β̂_n is obtained by maximizing the marginal likelihood dℙ(𝒟_n | β) with respect to β. Note that since f_0 and f_1 possess different smoothness parameters, then the “true" model for generating 𝒟_n is characterized by a likelihood function dℙ(𝒟_n | α_0, α_1). Assume that the true model dℙ(𝒟_n | α_0, α_1) is identifiable, i.e. the mapping (α_0, α_1) ↦ℙ is one-to-one. Type-I priors re-parametrize the observation model so that the likelihood function dℙ(𝒟_n | β) is parametrized with a single smoothness parameter β. Hence, as long as α_0 ≠α_1, the new parametrization renders an unidentifiable model, since the mapping β↦ℙ is not one-to-one (i.e. different combinations of α_0 and α_1 can map to the same β). This means that in this case likelihood-based empirical Bayes would never satisfy the matching condition in Theorem 4, even in the limit of infinite samples (n ↑∞). In most practical Bayesian models (e.g. Gaussian processes), the empirical Bayes estimate β̂_n will be in the interval (α_0,α_1) with high probability as depicted in Figure <ref>. This means that with high probability, the likelihood-based empirical Bayes method will prompt an oversmoothed prior, from which all draws are smoother than the true ITE function, leading to a suboptimal information rate. 

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 β̂_n that minimizes the posterior Bayesian KL divergence, i.e.   

    β̂_n    = min_β𝔼_Π(. | 𝒟_n,β)[ 𝔻_n(Π(β); f_0, f_1) | 𝒟_n ] 
       = min_β𝔼_Π(. | 𝒟_n,β)[ 𝔼_x[ D_ KL(P(x)  Q_𝒟_n(x)) ] ].
 
The information-based empirical Bayes estimator is simply a Bayesian estimator of β with the loss function being the posterior KL risk in (<ref>). Unlike the likelihood-based method, the objective in (<ref>) is an direct measure for the quality of causal inference conducted with a prior Π_β. In the following Theorem, we show that β̂_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^α_0 and H^α_1, respectively, and let Π(β) be a prior over Hölder space with order β. If β̂_n is obtained as in (<ref>) using cross-validation, then under certain regularity conditions we have that β̂_n p→ (α_0∧α_1).
  


Proof. See Appendix D.     


    
Theorem 5 says that the information-based empirical Bayes estimator is consistent. That is, the estimate β̂_n will eventually converge to α_0∧α_1 as n →∞. Note that this is a weaker result than adaptivity: consistency of β̂_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 α_0∧α_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 <cit.>. A more flexible Type-II prior that assigns different smoothness parameters β_0 and β_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 <cit.>. As a consequence of Theorem 1, it turns out that the information-based empirical Bayes estimator in (<ref>) is structurally similar to the risk-based empirical Bayes adaptation method proposed in <cit.>. Hence, we conjecture that our proposed empirical Bayes procedure can guarantee frequentist coverage for the estimated causal effects under some conditions <cit.>.




  §.§.§ 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 Π(β) = 𝒢𝒫((β)), which places a probability distribution over a Hölder space with regularity β, 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) <cit.>, 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))^-(α_0 ∧α_1)/2+d/4, whereas spline kernels achieved a rate of (n/log(n))^-2(α_0 ∧α_1)/(2(α_0 ∧α_1) + d), which is optimal up to a logarithmic factor <cit.>. 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 <cit.>. 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 (<ref>).   

Another possible option for Bayesian nonparametric priors place their probability mass on the space of piece-wise constant functions (trees) <cit.>. 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 Competition[http://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(α_0 ∧α_1)/(2α_0 ∧α_1 + d) <cit.>. 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 (<ref>).  

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.
 


§ 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 <ref>). Moreover, the RKHS space should be defined through a Matérn covariance kernel with parameters β_0 and β_1 for response surfaces f_0 and f_1 (Subsection <ref>), and the parameters β = (β_0,β_1) should be optimized via the information-based empirical Bayes procedure in (<ref>). In this Section, we construct a practical learning algorithm that follows this recipe.

Type-II GP priors place a probability distribution on functions f: 𝒳→ℝ^2 in a vector-valued Reproducing Kernel Hilbert Space (vvRKHS). A vvRKHS ℋ_ K is equipped with an inner product ⟨ .,. ⟩_ℋ_ K, and a reproducing kernel K: 𝒳×𝒳→ℝ^2 × 2, where K is a (symmetric) positive semi-definite matrix-valued function <cit.>. Note that by operating in a vvRKHS we get the algorithmic advantage of being able to conduct posterior inference in an infinite-dimensional function space by estimating a finite number of coefficients evaluated at the input feature points (this is a consequence of the well-known representer Theorem <cit.>). GP regression in vvRKHS is often associated with multi-task learning <cit.>, and the corresponding GP is often known as a multi-task GP <cit.>. Multi-task learning is a common setup in machine learning where one model shares parameters between different tasks to improve statistical efficiency. The results of Theorem 3 can be thought of as suggesting multi-task learning as a framework for causal inference, where learning each of the potential outcomes (f_0 and f_1) is thought of as a separate learning task, and a single model is used to execute the two tasks simultaneously. 

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 β_0 and β_1 to f_0 and f_1, respectively. Standard intrinsic coregionalization models for vector-valued kernels impose the same covariance parameters for all outputs <cit.>, 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) <cit.>, which mixes two intrinsic coregionalization models as follows 

    K_θ(x, x^') =  A k_0(x, x^') +  B  k_1(x, x^'),

where k_ω(x, x^') = (β_ω), ω∈{0,1}, whereas A and B are given by 

    A = [ a^2_00   a_01;   a_10      ϵ ],   B = [    ϵ b_01; b_10 b_11 ],

where ϵ→ 0 is a small positive number that is determined a priori and kept fixed during the prior adaptation procedure. The LCM kernel structure in (<ref>) ensures that the response surfaces f_0 and f_1 have smoothness levels β_0 abd β_1 respectively. The constant ϵ ensures that K_θ(x, x^') is positive semi-definite for any selection of the other parameters. The parameters a_00 and b_11 represent the variances of f_0 and f_1, whereas all other variables (a_01,a_10,b_01,b_10) are correlation variables that share information among the two learning tasks (learning f_0 and f_1). The set of all kernel parameters is denoted as β. Given a set of “hyper-parameters" β, the ITE function estimate T̂_n is obtained in terms of the posterior mean[Closed-form expressions for the posterior mean of a multi-task GP can be found in <cit.>.] as follows: T̂_n = 𝔼_Π_β[  f^T v | 𝒟_n], where v = [-1,1]^T.   

Now that we completely specified the multi-task GP prior for a given hyper-parameter set β, the only remaining ingredient in the recipe is to implement the information-based empirical Bayes adaptation criterion in (<ref>). The following Theorem gives an insightful decomposition of the information-based empirical Bayes objective for the multi-task GP model. (In the following Theorem, Y^(W) = [Y_i^(ω_i)]_i and Y^(1-W) = [Y_i^(1-ω_i)]_i are vectors comprising all factual and counterfactual outcomes associated with an observational dataset 𝒟_n.)



Theorem 6. (Factual bias and counterfactual variance decomposition) The minimizer β^* of the information-based empirical Bayes adaptation criterion in (<ref>) is given by

    min_β Y^(W)-𝔼_Π_β[  f | 𝒟_n ]_2^2_ + _Π_β[  Y^(1-W) | 𝒟_n ]_1_,

where _Π_β 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 (<ref>) decomposes to factual bias and counterfactual variance terms[The objective function in Theorem 6 can be easily optimized via a leave-one-out cross-validation procedure. Refer to <cit.> 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 β. In that sense, the factual bias is a measure of the goodness-of-fit for the posterior mean resulting from a prior smoothness β. 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 β. A small value for β 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 β 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 α_0 and α_1, and estimating the value of  α_0 ∧α_1. That is, the regularizer attempts to protect the prior from falsely recognizing either α_0 or α_1 as being very low just because it over-fit the factual outcomes, and hence underestimating the true α_0 ∧α_1, thereby undersmoothing the prior and giving rise to a suboptimal information rate. The two terms work in opposite directions as shown in Figure <ref>: 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.   






  


   











§ EXPERIMENTS


We sought to evaluate the finite-sample performance of the Bayesian causal inference procedure proposed in Section <ref>, and compare it with state-of-the-art causal inference models. Causal inference models are hard to evaluate <cit.>, 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 <cit.>. 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 IHDP dataset

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 <cit.>. The semi-simulated dataset in <cit.> 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 <cit.>. The dataset comprises 747 subjects (608 control and 139 treated), and there are 25 features associated with each subject.

 



 §.§ Benchmarks
 
We compared our algorithm with various causal models and standard machine learning benchmarks which we list in what follows:  Tree-based methods (BART <cit.>, causal forests (CF) <cit.>,  Balancing counterfactual regression (balancing neural networks (BNN) <cit.>, and counterfactual regression with Wasserstein distance metric (CFRW) <cit.>),  Propensity-based and matching methods (k nearest-neighbor (kNN), propensity score matching (PSM)), a  nonparametric spline regression model (causal MARS <cit.>), and ⊙ Doubly-robust methods (Targeted maximum likelihood (TML) <cit.>). 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 □ Type-I regression structure, in which we use the treatment assignment variable as an input feature to the machine leaning algorithm, and a ⊗ 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 <cit.>, and did not tune the model's hyper-parameters. For BNN and CFRW, we used the neural network configurations reported in <cit.> and <cit.>. Causal MARS was implemented as described in <cit.>. PSM was implemented as described in <cit.>, 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 , with the exception of BART, causal forests and TMLE, all of which were implemented in . We used the  libraries , , and  for the implementation of BART, causal forests and TMLE, respectively. Our method was implemented in  using , a library for Gaussian processes <cit.>.      



 §.§ Evaluation
 
We evaluate the performance of all benchmarks by reporting the square-root of the PEHE. The empirical PEHE is estimated as = 1/n∑_i=1^n((f_1(X_i)-f_0(X_i))-𝔼[Y^(1)_i-Y^(0)_i|X_i=x])^2, where f_1(X_i)-f_0(X_i) is the estimated treatment effect. We evaluate the PEHE of all algorithms via a Monte Carlo simulation with 1000 realizations of the IHDP semi-synthetic model, where in each experiment/realization we run all the benchmarks with a 60/20/20 train-validation-test splits. (For models that do not need hyper-parameter tuning, such as BART and our GP models, the entire training set is used to compute the posterior distributions.) We report both the in-sample and out-of-sample PEHE estimates: the former corresponds to the accuracy of the estimated ITE in a retrospective cohort study, whereas the latter corresponds to the performance of a clinical decision support system that provides out-of-sample patients with ITE estimates <cit.>. The in-sample PEHE results are non-trivial since we never observe counterfactuals even in the training phase. Recall that, from Theorem 1, we know that the achieved information rate by a Bayesian inference procedure is equivalent to the PEHE estimation rate. Thus, the PEHE performance is a direct proxy of the achieved information rate, and since it is an essentially frequentist quantity, we can use it to compare the performance of our model with the frequentist benchmarks.         



 §.§ Results
 
As can be seen in Table <ref>, 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 <ref> 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 <cit.>). Our algorithm is also performing better than all other nonparametric tree-based algorithms. This is expected since, as we have discussed earlier in Subsection <ref>, 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 <cit.>), and hence all machine learning methods exhibit high generalization errors. 



§ 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. 




§ PROOF OF THEOREM 2


We start by establishing an asymptotic equivalence between the KL risk and the frequentist loss in the L_2(ℙ) norm, i.e. 

    𝔻_n(Π; f_0, f_1,γ) ≍𝔼_𝒟_n[ ‖𝔼_Π[ T | 𝒟_n ]-T‖^2_2 ].

Note that, from Theorem 1, we already know that since the expected KL risk is bounded above by the L_2(ℙ) loss (with a constant factor), then it follows that: 

    𝔻_n(Π; f_0, f_1,γ) ≲𝔼_𝒟_n[ ‖𝔼_Π[ T | 𝒟_n ]-T‖^2_2 ].

Recall from (<ref>) that the KL risk is given by

    𝔻_n(Π; f_0, f_1,γ) = 𝔼_x[ 𝔼_𝒟_n[D_ KL(P(x)  Q_𝒟_n(x)).] ].

Using Pinsker's inequality <cit.>, the KL divergence can be bounded below as follows:

    P(x)-Q_𝒟_n(x)_TV≤√(1/2 D_ KL(P(x)  Q_𝒟_n(x)).),

where ._TVis the total variation distance between probability measures, which is given by the L_1 norm of the difference between P(x) and Q_𝒟_n(x) as follows:

    P(x)-Q_𝒟_n(x)_TV   = P(x)-Q_𝒟_n(x)_1
       = ∫_𝒳 |P(x)-Q_𝒟_n(x)|  dx.
 
Since the L_1 norm is bounded below by the L_2 norm, we can lower bound the KL divergence by combining (<ref>) with Pinsker's inequality as follows: 

    D_ KL(P(x)  Q_𝒟_n(x)).    ≥ 2 P(x)-Q_𝒟_n(x)^2_1
       ≥ 2 P(x)-Q_𝒟_n(x)^2_2,
 
and hence it follows that

    𝔻_n(Π; f_0, f_1,γ)    ≥ 2 𝔼_x[ 𝔼_𝒟_n[P(x)-Q_𝒟_n(x)^2_2] ] 
       = 2 𝔼_x[ 𝔼_𝒟_n[T-𝔼_Π[ T | 𝒟_n ]^2_2] ],

which leads to the following asymptotic inequality

    𝔻_n(Π; f_0, f_1,γ) ≳𝔼_𝒟_n[T-𝔼_Π[ T | 𝒟_n ]^2_2].

By combining (<ref>) and (<ref>), we arrive at (<ref>). From (<ref>), it follows that the optimal information rate is equivalent to the minimax estimation rate in the L_2(ℙ) norm, i.e.

    I^*_n(ℱ_0, ℱ_1,γ) ≍min_T̂max_f_0 ∈ℱ_0, f_1 ∈ℱ_1𝔼_𝒟_n[T-T̂^2_2],

where the estimator T̂ is taken to be 𝔼_Π[ T | 𝒟_n ] since the posterior mean estimator is optimal for the L_2(ℙ) loss. In what follows, we derive the optimal information rate by obtaining the minimax rate of estimation the ITE function T(x).     

Let δ_ω be the solution to H(δ_ω; ℱ^α_ω) ≍ n δ^2_ω. We will prove that the optimal rate is Θ(δ^2_0 ∨δ^2_1) by first showing that I^*_n(ℱ^α_0,ℱ^α_1) is lower bounded by, i.e. I^*_n(ℱ^α_0,ℱ^α_1) = Ω(δ^2_0 ∨δ^2_1), and then show that I^*_n(ℱ^α_0,ℱ^α_1) = O(δ^2_0 ∨δ^2_1). We start by observing that the causal inference problem can be described through the following Markov chain 

    (f_0,f_1) →𝒟_n → (f̂_0,f̂_1) →T̂.

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

    I(T;T̂)  (*)≤ I(T;𝒟_n)  (⋆)≤sup_Π I(T;𝒟_n),

where (*) follows from the data processing inequality <cit.>, and the supremum in (⋆) is taken over all possible priors. I(T;T̂) is bounded below by the rate-distortion function

    I(T;T̂) ≥inf_T, T̂: 𝔼T-T̂_2^2  ≤  R_Π^* I(T;T̂),

for any T̂ satisfying 𝔼T-T̂_2^2 ≤ R_Π^*, where the infimum is taken over all joint distributions of (T,T̂). Combining (<ref>) and (<ref>), we can upper and lower bound the mutual information I(T;T̂) as follows

    inf_ 𝔼T-T̂_2^2  ≤  R_Π^* I(T;T̂) ≤ I(T;T̂) ≤sup_Π I(T;𝒟_n).

The lower bound in the chain of inequalities above is intractable, and hence we further lower bound I(T;T̂) using Fano's method <cit.>. That is, we take discrete subsets ℱ̃^α_0 and ℱ̃^α_1 of the function spaces ℱ^α_0 and ℱ^α_1, and convert the estimation problem to a testing problem. The spaces 

    ℱ̃^α_ω = {f̃^1_ω,. . .,f̃^M̃_ω_ω},  ℱ̃^α_ω⊂ℱ^α_ω,  ω∈{0,1},

are constructed such that f̃^i_ω-f̃^j_ω≥δ, ∀ i ≠ j. Let Q be a quantizer that maps elements of ℱ^α_ω to ℱ̃^α_ω,  ω∈{0,1}. Thus, the causal inference problem can be described through the following Markov chain: 

    (f_0,f_1) →𝒟_n → (f̂_0,f̂_1) → Q(f̂_0,f̂_1).

Let T̃ = f̃^u_1-f̃^v_0, where f̃^v_0 and f̃^u_1 are the functions in ℱ̃^α_0 and ℱ̃^α_1 that are closest to f_0 and f_1. The discrete element T̃ belongs to a set {T̃^1,. . .,T̃^M̃_T}, which corresponds to a discretized version of the function space to which T belongs. Using the data processing inequality, we have that  

    I(T̃;T̂) ≥ I(T̃;Q(T̂)).

An “error event" is an event where Q(T̂) does not correspond to the true discretized function T̃, i.e. the event {T̃≠ Q(T̂)}. The error event occurs when

    T̂-Q(T̂)≤T̂-T̃, {T̃≠ Q(T̂)}.

Thus, the error event implies that δ≤Q(T̂)-T̃. Using the triangular inequality, (<ref>) can be further bounded as follows: 

    δ   ≤Q(T̂)-T̃ = Q(T̂)-T̂+T̂-T̃
       ≤Q(T̂)-T̂ + T̂-T̃
       ≤ 2 T̂-T̃ T̂-T̃≥δ/2.

Let P_e be the probability of the error event {T̃≠ Q(T̂)}. From (<ref>), P_e can be bounded above as follows 

    P_e    ℙ({T̃≠ Q(T̂)}) 
       = ℙ(Q(T̂)-T̃≥δ) = ℙ(T̂-T̃≥δ/2) 
       = ℙ(T̂-T̃_2^2 ≥δ^2/4) 
        (∙)≤4/δ^2 𝔼[T̂-T̃_2^2] ≤4/δ^2 R^*_Π,

where (∙) is an application of Markov's inequality. By combining (<ref>) with the result in (<ref>), the lower bound in (<ref>) can be further bounded below as follows 

    inf_𝔼T-T̂_2^2  ≤  R_Π^* I(T;T̂)    ≥inf_𝔼T-T̂_2^2  ≤  R_Π^* I(T̃;T̂) 
       = inf_P_e ≤4/δ^2 R^*_Π I(T̃;T̂) 
       ≥inf_P_e ≤4/δ^2 R^*_Π I(T̃;Q(T̂)).

The mutual information I(T̃;Q(T̂)) can be bounded above as follows

    I(T̃;Q(T̂))    = I(f̃_1-f̃_0;Q(f̂_1-f̂_0))
        (⊙)≤ I(f̃_0,f̃_1;Q(f̂_1-f̂_0))
       ≤ I(f̃_0,f̃_1;Q(f̂_0),Q(f̂_1))
       = I(f̃_0;Q(f̂_0)) + I(f̃_1; Q(f̂_1))
       ≤ 2 {I(f̃_0;Q(f̂_0)), I(f̃_1;Q(f̂_1))},

where (⊙) follows from the data processing inequality. Note that the mutual information I(T̃;Q(T̂)) can be written in terms of the KL divergence as <cit.> 

    I(T̃;Q(T̂))    = D(ℙ(T̃;Q(T̂)) || ℙ(T̃) ·ℙ(Q(T̂))) 
       ≥ D((P_e) || (1-1/n)) 
       = P_elog(P_e/1-1/M̃_T) + (1-P_e)log(1-P_e/1/M̃_T) 
       = -h(P_e) + log(M̃_T)-P_e log(M̃_T-1) 
       ≥ -log(2) + log(M̃_T) - P_e log(M̃_T),

where h(.) is the binary entropy. From (<ref>), we have that 

    P_e    ≥ 1-I(T̃;Q(T̂)) + log(2)/log(M̃_T),

which is an incarnation of Fano's inequality. By combining (<ref>) with (<ref>), we have the following inequality

    P_e    ≥ 1-I(f̃_0;Q(f̂_0)) ∨ I(f̃_1;Q(f̂_1)) + log(√(2))/1/2 log(M̃_T).

From (<ref>), the minimax risk R^*_Π is bounded below by

    R^*_Π   ≥δ^2/4(1-I(f̃_0;Q(f̂_0)) ∨ I(f̃_1;Q(f̂_1)) + log(√(2))/1/2 log(M̃_T)).

The discretization ℱ̃^α_ω = {f̃^1_ω,. . .,f̃^M̃_ω_ω} corresponds to a δ-packing of the function space ℱ^α_ω, and hence M̃_ω is given by the covering number N(δ, ℱ^α_ω), for ω∈{0,1}. It follows that M̃_T ≥ N(δ, ℱ^α_0) ∨ N(δ, ℱ^α_1), and hence we have that

    R^*_Π   ≥δ^2/4(1-I(f̃_0;Q(f̂_0)) ∨ I(f̃_1;Q(f̂_1)) + log(√(2))/1/2 log(N(δ, ℱ^α_0) ∨ N(δ, ℱ^α_1))).

The mutual information I(f̃_ω;Q(f̂_ω)) can be bounded via the KL divergence as

    I(f̃_ω;Q(f̂_ω))    ≤1/N^2(δ, ℱ^α_ω)∑_i,j D(ℙ(f̃^i_ω) || ℙ(f̃^j_ω)) 
       ≤ 2 n δ^2.

Thus, the minimax risk can be bounded below as follows

    R^*_Π   ≥δ^2/4(1-4 n δ^2 + log(2)/log(N(δ, ℱ^α_0) ∨ N(δ, ℱ^α_1))),

and hence we have that 

    R^*_Π   ≳δ^2 - δ^4 n + δ^2/log(N(δ, ℱ^α_0) ∨ N(δ, ℱ^α_1)).

Since R^*_Π is strictly positive, then we have that 

    R^*_Π   ≳δ^2,

where δ is the solution to the transcendental equation

    δ^2 ≍δ^4 n/log(N(δ, ℱ^α_0) ∨ N(δ, ℱ^α_1)),

or equivalently

    log(N(δ, ℱ^α_0) ∨ N(δ, ℱ^α_1)) ≍δ^2 n.

The metric entropy of a function space ℱ^α_ω is given by H(δ, ℱ^α_ω) = log(N(δ, ℱ^α_ω), and hence (<ref>) is written as 

    H(δ, ℱ^α_0) ∨ H(δ, ℱ^α_1) ≍δ^2 n.

Since the metric entropy H(δ, ℱ^α_ω) is a decreasing function of the smoothness parameter α_ω, then it follows that the solution δ^* of the transcendental equation in (<ref>) is given by δ^* = δ_0 ∨δ_1, where δ_ω is the solution to the equation   

    H(δ_ω, ℱ^α_ω) ≍δ_ω^2 n, ω∈{0,1}.

The equation in (<ref>) has a solution for all n when the function space ℱ^α_ω has a polynomial or a logarithmic metric entropy <cit.>, which is the case for all function spaces of interest (see Table I for evaluations of δ_0 ∨δ_1 for various function spaces). It follows from (<ref>) and (<ref>) that

    R^*_Π = Ω(δ^2_0 ∨δ^2_1),   H(δ_ω, ℱ^α_ω) ≍δ_ω^2 n,  ω∈{0,1},

and hence, from (<ref>), we have that

    I^*_n = Ω(δ^2_0 ∨δ^2_1),   H(δ_ω, ℱ^α_ω) ≍δ_ω^2 n,  ω∈{0,1}.

We now focus on upper bounding R^*_Π. From <cit.>, we know that the minimax risk is upper bounded by the channel capacity in (<ref>), which is further bounded above by the covering numbers as follows 

    R^*_Π≲1/n( log(N(δ,ℱ^α_0)) ∨log(N(δ,ℱ^α_1)) + n δ^2 ).

For δ satisfying (<ref>), we have that

    log(N(δ,ℱ^α_0)) ∨log(N(δ,ℱ^α_1)) = δ^2 n,

and hence R^*_Π≲δ^2_0 ∨δ^2_1. It follows that        

    I^*_n = O(δ^2_0 ∨δ^2_1),   H(δ_ω, ℱ^α_ω) ≍δ_ω^2 n,  ω∈{0,1}.

By combining (<ref>) and (<ref>), we have that I^*_n = Ω(δ^2_0 ∧δ^2_1) and I^*_n = O(δ^2_0 ∨δ^2_1), and hence it follows that 

    I^*_n = Θ(δ^2_0 ∨δ^2_1),   H(δ_ω, ℱ^α_ω) ≍δ_ω^2 n,  ω∈{0,1}.




§ PROOF OF THEOREM 3


Note that when f_0 ∈ H^α_0_𝒫_0 and f_1 ∈ H^α_1_𝒫_1, the metric entropy of H^α_0_𝒫_0 and H^α_1_𝒫_1 are given by <cit.>: 

    H(δ,H^α_0_𝒫_0) ≍δ^-|𝒫_0|/α_0,   H(δ,H^α_1_𝒫_1) ≍δ^-|𝒫_1|/α_1.

From Theorem 2, we know that the optimal information rate is given by I^*_n(H^α_0_𝒫_0,H^α_1_𝒫_1) ≍δ_0^2 ∨δ_1^2, where δ_0 and δ_1 are the solutions for δ^-|𝒫_ω|/α_ω≍ n δ^2_ω, ω∈{0,1}. Thus, we have that

    δ_ω≍ n^-2α_ω/2α_ω + |𝒫_ω|, ω∈{0,1},

and hence the optimal information rate is given by  

    I^*_n(H^α_0_𝒫_0,H^α_1_𝒫_1) ≍ n^-2α_0/2α_0+|𝒫_0|∨ n^-2α_1/2α_1+|𝒫_1|,

From Theorem 1, we know that

    𝔻_n(Π; f_0, f_1,γ)    ≍𝔼[ ‖ 𝔼_Π[ T | 𝒟_n ]-T ‖_2^2 ].

The L_2(ℙ) loss term on the right hand side can be upper bounded as follows (see Lemma C2): 

    𝔼[ ‖ 𝔼[ T | 𝒟 ]-T ‖_2^2 ]    ≲𝔼[ ‖ 𝔼[ f_0 | 𝒟 ]-f_0 ‖_2^2 ] 
       +  𝔼[ ‖ 𝔼[ f_1 | 𝒟 ]-f_1 ‖_2^2 ],

For a Type-II prior Π∈Π_β_0,1^∘∘ over the two Hölder spaces H^α_0_𝒫_0 and H^α_1_𝒫_1, with β_0 = α_0 and β_1 = α_1, the minimax estimation rates for nonparametric regression over f_0 and f_1 are <cit.>

    inf_f̂_0sup_f_0 ∈ H^α_0𝔼[  f̂_0(𝒟_n)-f_0 ‖_2^2 ]    ≍ n^-2α_0/2α_0+|𝒫_0|,
    inf_f̂_1sup_f_1 ∈ H^α_1𝔼[  f̂_1(𝒟_n)-f_1 ‖_2^2 ]     ≍ n^-2α_1/2α_1+|𝒫_1|,

and it follows that

    ℐ_n(Π; f_0, f_1,γ) = O(n^-2α_0/2α_0+|𝒫_0|∨ n^-2α_1/2α_1+|𝒫_1|),

which matches the optimal information rate in Theorem 2. Similarly, for a Type-I prior Π∈Π_β^∘ over a Hölder space H^β_𝒫_0 ∪𝒫_1, with β = α_0 ∧α_1, the term 𝔼[ ‖ 𝔼[ f_ω | 𝒟 ]-f_ω ‖_2^2 ] becomes equivalent to the L_2(ℙ) of nonparametric regression of the surface f : [0,1]^|𝒫_0 ∪𝒫_1|×{0,1}→ℝ. The minimax estimation rate of such a problem is <cit.>

    inf_f̂sup_f ∈ H^α_0𝔼[  f̂(𝒟_n)-f ‖_2^2 ]    ≍ n^-2(α_0 ∧α_1)/2β+|𝒫_0 ∪𝒫_1|,
 
where the number of feature dimensions |𝒫_0 ∪𝒫_1| correspond to all the relevant dimensions for the regression function f(x,ω). The regression function on the discrete dimension ω can be estimated at the √(n) parametric rate and hence it does not affect the minimax estimation rate given above. Since the rate n^-2(α_0 ∧α_1)/2β+|𝒫_0 ∪𝒫_1| is strictly slower than the optimal rate of n^-2α_0/2α_0+|𝒫_0|∨ n^-2α_1/2α_1+|𝒫_1| for all β > 0, it follows that 

    𝔻_n(Π; f_0, f_1,γ) ≲ I^*_n(H^α_0_𝒫_0,H^α_1_𝒫_1), ∀Π∈Π_β^∘, ∀β > 0.
 



§ PROOF OF THEOREM 4


We start by providing three Lemmas which we will use to prove the statement of the Theorem.



Lemma C1. Let 𝒳 be a compact subset of ℝ^d, α, β∈ [0,1], and n, m ∈ℕ_0. If n + β > m + α, then H^n+β(𝒳) is compactly contained in H^m+α(𝒳).     



Lemma C2. The L_2(ℙ) loss L = 𝔼[ ‖ 𝔼[ T | 𝒟_n ]-T ‖_2^2 ] is asymptotically bounded above as follows:

    L    ≲𝔼[ ‖ 𝔼[ f_0 | 𝒟_n ]-f_0 ‖_2^2 ] + 𝔼[ ‖ 𝔼[ f_1 | 𝒟_n ]-f_1 ‖_2^2 ].
 
Proof. The L_2(ℙ) loss conditioned on an observational dataset L(𝒟_n) = ‖ 𝔼[ T | 𝒟_n ]-T ‖_2^2 is given by:

    L(𝒟_n)=  ‖ (f̂_1(x)-f̂_0(x))-(f_1(x)-f_0(x)) ‖_2^2,
 
where f̂_ω(x) = 𝔼[ f_ω(x) | 𝒟_n ], for ω∈{0,1}. The L_2(ℙ) norm in (<ref>) can be expressed as follows:

    L(𝒟_n)    =  ‖ (f̂_1(x)-f̂_0(x))-(f_1(x)-f_0(x)) ‖_2^2, 
       = ∫_𝒳 ((f̂_1(x)-f̂_0(x))-(f_1(x)-f_0(x)))^2  dℙ(x) 
       = ∫_𝒳 ((f̂_1(x)-f_1(x))+(f_0(x)-f̂_0(x)))^2  dℙ(x) 
       ≤ 2 ∫_𝒳 ((f̂_1(x)-f_1(x))^2+(f̂_0(x)-f_0(x))^2)  dℙ(x) 
       = 2 ∫_𝒳 (f̂_1(x)-f_1(x))^2  dℙ(x,ω=1) 
       + 2 ∫_𝒳 (f̂_0(x)-f_0(x))^2  dℙ(x,ω=0).
 
Since dℙ(x,ω=1) = γ(x)· dℙ(x) and dℙ(x,ω=0) = (1-γ(x))· dℙ(x), we have that

    L(𝒟_n) 
       = 2 ∫_𝒳 (f̂_1(x)-f_1(x))^2 γ(x)· dℙ(x) 
       + 2 ∫_𝒳 (f̂_0(x)-f_0(x))^2  (1-γ(x))· dℙ(x)
       = 2 √(γ(x))·(f̂_1(x)-f_1(x))^2_L_2(ℙ)
       + 2 √(1-γ(x))·(f̂_0(x)-f_0(x))^2_L_2(ℙ).

Using Cauchy-Schwarz inequality, we obtain the following:

    √(γ(x))(f̂_1(x)-f_1(x))^2_2 ≤γ(x)_2 ·(f̂_1(x)-f_1(x))^2_2,

and similarly for √(1-γ(x))·(f̂_0(x)-f_0(x))^2_2. The proof of the Lemma is concluded by observing that γ(x)_2 is O(1) and (f̂_1(x)-f_1(x))^2_2 ≍(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 Π(β) = 𝒢𝒫((β)) is the space of Hölder functions with order β.    



The proofs for Lemmas C1 and C3 are standard and can be found in <cit.> and <cit.> respectively.
  
Recall that the expected Kullback-Leibler risk and the L_2(ℙ) loss are asymptotically equivalent (see Appendix A):

    𝔻_n(Π; f_0, f_1,γ) ≍𝔼_𝒟_n[ ‖𝔼_Π[ T | 𝒟_n ]-T‖^2_2 ].

From Lemma C2 and the equivalence in (<ref>), we have that

    I_n(Π(β_0,β_1); H^α_0, H^α_1)    ≍𝔼_𝒟_n[ ‖𝔼_Π[ T | 𝒟_n ]-T‖^2_2 ]. 
       ≲𝔼_𝒟_n[ ‖ 𝔼_Π[ f_0 | 𝒟 ]-f_0 ‖_2^2 ] 
       + 𝔼_𝒟_n[ ‖ 𝔼_Π[ f_1 | 𝒟 ]-f_1 ‖_2^2 ].

Thus, the information rate achieved by the prior Π(β_0,β_1) is upper bounded by the posterior contraction rates <cit.> (rate of convergence of the L_2(ℙ) loss) over the surfaces f_0 and f_1. For a prior 𝒢𝒫((β_ω)) and a true function f_ω∈ H^α_ω, the contraction rate ε^2 is given by solving the following transcendental equation  <cit.>:    

    ϕ_f_ω(ε) ≍ n·ε^2,

where ϕ_f_ω(ε) is the concentration function defined as <cit.>:

    ϕ_f_ω(ε)  - log(ℙ_Π(β_ω)(f-f_ω_∞ < ε)).
 
The concentration function measures the amount of prior mass that Π places around the true function f_ω. The transcendental equation in (<ref>) provides a valid contraction rate whenever consistency holds. Consistency of Bayesian inference holds whenever the true parameter (in this case the true function f_ω) is in the support of the prior <cit.>. From Lemmas C1 and C3, it follows that the necessary and sufficient conditions for consistency is that β_0 ≤α_0 and β_1 ≤α_1.  

In <cit.>, the concentration function ϕ_f_ω(ε) for a sufficiently smooth prior 𝒢𝒫((β_ω)), with β_ω > d/2, and a sufficiently smooth true function f_ω∈ H^α_ω, with β_ω > d/2, was obtained as follows: 

    ϕ_f_ω(ε) ≲ε^-d/β_ω + ε^-2β_ω-2α_ω +d/α_ω.
  
Thus, combining (<ref>) and (<ref>), the posterior contraction rate for Π(β_ω) around f_ω is the solution to:  

    n·ε^2 ≲ε^-d/β_ω + ε^-2β_ω-2α_ω +d/α_ω,

The solution to (<ref>) is given by

    ε   ≲ n^-β_ω/2β_ω + d + n^-α_ω/2β_ω+d, 
       ≍ n^-(β_ω∧α_ω)/2β_ω + d.

Since consistency holds for β_ω≤α_ω, then ε = O(n^-β_ω/2β_ω + d) and the contraction rate is n^-2β_ω/2β_ω + d. That is, we can characterize the L_2(ℙ) loss surfaces on f_0 and f_1 as follows:

    𝔼_𝒟_n[ ‖ 𝔼_Π[ f_ω | 𝒟 ]-f_ω ‖_2^2 ] ≲ n^-2β_ω/2β_ω + d,

and so it follows that:

    I_n(Π(β_0,β_1); H^α_0, H^α_1) ≲ n^-2β_0/2β_0 + d + n^-2β_1/2β_1 + d ,

which concludes the proof of the Theorem. 



§ PROOF OF THEOREM 5


The empirical smoothness estimate β̂_n is obtained by minimizing the empirical objective:

    L(𝒟_n, β) = 𝔼_f_0,f_1 ∼ dΠ_β(. | 𝒟_n)[ 𝔻_n(Π_β; f_0, f_1) | 𝒟_n ].

We optimize (<ref>) via model selection with J-fold cross-validation. From Theorem 4, we know that for a large sample (n ↑∞), the true loss function L(β) = 𝔼_𝒟_n[L(𝒟_n, β)] has a unique minimizer: β^* = (α_0 ∧α_1). Let the set {β^(1),. . .,β^(K_n)} be a set of K_n candidate minimizers (smoothness levels) of the true loss L. Let B_n = {B_n(i)}_i ∈{0,1}^n be a binary split vector which allocates every data point i in 𝒟_n to either of the training or validation sets. We define ℙ^T_n,B_n and ℙ^V_n,B_n be the empirical distributions of the training and validation sets, and let v be the fraction of data allocated to the validation set. The empirical cross-validated risk estimate is defined as

    L^k_n(𝒟_n, β^(k_n))    = 𝔼_B_n∫ L(𝒟_n,ℙ^T_n,B_n, β^(k_n))  dℙ^V_n,B_n
       = 𝔼_B_n  L(𝒟^i_n,ℙ^T_n,B_n).

The candidate in {β^(1),. . .,β^(K_n)} that minimizes the cross-validated risk L^k_n(𝒟_n) is 

    k̂^*_nmin_k_n L(𝒟_n,β^(k_n)).

The consistency of the estimator β^(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 <cit.>. Suppose that sup_𝒟_n,β L(𝒟_n,β) ≤∞, β^* ∈{β^(1),. . .,β^(K_n)}, and log(K_n)/(√(nv)(L^k_n-L)) p→ 0 as n →∞. Then, from Theorem 2 in <cit.>, we have that L^k_n(𝒟_n, β^(k_n))-L(𝒟_n,β) p→ 0. Since β^* = (α_0 ∧α_1) is a unique minimizer of L(𝒟_n,β), then it follows from the argmin continuous mapping theorem for M-estimators that β^(K^*_n)→ (α_0 ∧α_1) <cit.>.   



§ PROOF OF THEOREM 6


Using the asymptotic equivalence between the KL and the L_2(ℙ) risks, we have that 

    β̂_n    = min_β𝔼_Π(. | 𝒟_n,β)[ 𝔼_x[ D_ KL(P(x)  Q_𝒟_n(x)) ] ] 
       ≍min_β𝔼_Π(. | 𝒟_n,β)[ 𝔼_x[ 𝔼_Π[T | 𝒟_n]-T^2_2 ] ].
 
Since there is a unique minimizer for both the KL and L_2(ℙ) risks when the number of samples is asymptotically large, we work on the L_2(ℙ) risk objective to obtain the optimal solution. The posterior Bayesian risk R(θ,f̂; 𝒟) for a point estimate f̂ is given by        

    R(θ,f̂; 𝒟) = 𝔼_θ[.ℒ̂(f̂; K_θ, Y^(W), Y^(1-W)) | 𝒟],

where the expectation in is taken with respect to Y^(1-W)|𝒟. The Bayesian risk can be written as

    R(θ,f̂; 𝒟) = ∫ℒ̂(f̂; K_θ, Y^(W), Y^(1-W)) dℙ_θ( Y^(1-W)|𝒟).

The loss function ℒ̂ conditional on a realization of the counterfactual outcomes is given by 

    ℒ̂(f̂; K_θ, Y^(W), Y^(1-W)) = 
       1/n∑_i=1^n(f̂^T(X_i) e-(1-2W_i)(Y^(1-W_i)_i-Y^(W_i)_i))^2.



The optimal hyper-parameter and interpolant (f̂^*,θ^*) are obtained through the following optimization problem in (F.1). The optimization problem can solved separately for θ and f̂; we know from Theorem 1 that for any given θ, the optimal interpolant f̂ = 𝔼_θ[ f | 𝒟]. Hence, the optimal hyper-parameter θ^* can be found by solving the optimization problem in (F.2). The objective function R can thus be written as in (F.3) and further reduced as in (F.4).
 
Note that since Y_i^(W_i) = f_W_i(X_i) + ϵ_i,W_i, then we have that 𝔼_θ[f_W_i(X_i) | 𝒟] = 𝔼_θ[Y_i^(W_i) | 𝒟] and 𝔼_θ[f_1-W_i(X_i) | 𝒟] = 𝔼_θ[Y_i^(1-W_i) | 𝒟]. Therefore, we can evaluate the terms R_1, R_2 and R_3 as follows

    R_1    = 1/n∑_i=1^n∫ (Y^(W_i)_i-𝔼_θ[f_W_i(X_i) | 𝒟])^2   dℙ_θ(Y_i^(1-W)|𝒟) 
       = 1/n∑_i=1^n∫ (Y^(W_i)_i-𝔼_θ[Y^(W_i)_i | 𝒟])^2   dℙ_θ(Y_i^(1-W)|𝒟) 
       = 1/n Y^(W)-𝔼_θ[ f | 𝒟]_2^2,
 
and

    R_2    = 1/n∑_i=1^n∫ (Y^(1-W_i)_i-𝔼_θ[f_1-W_i | 𝒟])^2  dℙ_θ(Y_i^(1-W)|𝒟) 
       = 1/n∑_i=1^n∫ (Y^(1-W_i)_i-𝔼_θ[Y^(1-W_i)_i | 𝒟])^2  dℙ_θ(Y_i^(1-W)|𝒟) 
       = 1/n∑_i=1^n[  Y^(1-W_i)_i | 𝒟 ], 
       = 1/n[  Y^(1-W) | 𝒟 ]_1,
 
and

    R_3 =    1/n∑_i=1^n∫ (Y^(W_i)_i-𝔼_θ[f_W_i | 𝒟]) 
       (Y^(1-W_i)_i-𝔼_θ[f_1-W_i | 𝒟])   dℙ_θ(Y_i^(1-W)|𝒟) = 0
 
Therefore, θ^* is found by minimizing Y^(W)-𝔼_θ[ f | 𝒟]_2^2+[  Y^(1-W) | 𝒟 ]_1.


IEEEtran



