Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
68,712 characters · 18 sections · 85 citation commands
Double Robust Bayesian Inference on Average Treatment Effects
{ Keywords: Average treatment effects, unconfoundedness, double robustness, nonparametric Bayesian inference, Bernstein–von Mises theorem, Gaussian processes. }
This paper proposes a double robust Bayesian approach for estimating the average treatment effect (ATE) under unconfoundedness, given a set of pretreatment covariates. Our new Bayesian procedure involves both prior and posterior adjustments. First, following ray2020causal, we adjust the prior distributions of the conditional mean function using an estimator of the propensity score. Second, we use this propensity score estimator together with a pilot estimator of the conditional mean to correct the posterior distribution of the ATE. The adjustments in both steps are closely related to the functional form of the semiparametric influence function for ATE estimation under unconfoundedness. They do not only shift the center but also change the shape of the posterior distribution. For our robust Bayesian procedure, we derive a new Bernstein–von Mises (BvM) theorem, which means that this posterior distribution, when centered at any efficient estimator, is asymptotically normal with the efficient variance in the semiparametric sense. The key innovation of our paper is that this result holds under double robust smoothness assumptions within the Bayesian framework.
Despite the recent success of Bayesian methods, the literature on ATE estimation is predominantly frequentist-based. For the missing data problem specifically, it was shown that conventional Bayesian approaches (i.e., using uncorrected priors) can produce inconsistent estimates, unless some unnecessarily strong smoothness conditions on the underlying functions were imposed; see the results and discussion in robins1997 or ritov2014. Once the prior distribution was adjusted using some pre-estimated propensity score, ray2020causal recently established a novel semiparametric BvM theorem under weaker smoothness requirement for the propensity score function.\footnote{Strictly speaking, the main objective in ray2020causal concerns the mean response in a missing data model, which is equivalent to observing one arm (either the treatment or control) of the causal setup. } However, a minimum differentiability of order $p/2$ is still required for the conditional mean function in the outcome equation, where $p$ denotes the dimensionality of covariates. In this paper, we are interested in Bayesian inference under double robustness that allows for a trade-off between the required levels of smoothness in the propensity score and the conditional mean functions.
Under double robust smoothness conditions, we show that Bayesian methods, which use propensity score adjusted priors as in ray2020causal, satisfy the BvM Theorem only up to a “bias term” depending on the unknown true conditional mean and propensity score functions. In this paper, our robust Bayesian approach accounts for this bias term in the BvM Theorem by considering an explicit posterior correction. Both the prior adjustment and the posterior correction are based on functional forms that are closely related to the efficient influence function for the ATE, see hahn1998role. We show that the corrected posterior satisfies the BvM Theorem under double robust smoothness assumptions. Our novel procedure combines the advantages of Bayesian methodology with the robustness features that are the strengths of frequentist procedures. Our credible intervals are Bayesianly justifiable in the sense of rubin1984bayesian, as the uncertainty quantification is made conditional on the observed data and can be also interpreted as frequentist confidence intervals with asymptotically exact coverage probability. Our procedure is inspired by insights from the double machine learning (DML) literature, as well as the bias-corrected matching approach from abadie2011bias, as our robustification of an initial procedure removes some non-negligible bias and remains asymptotically valid under weaker regularity conditions. While the main part of our theoretical analysis focuses on the ATE of binary outcomes, also considered by ray2020causal, we outline extensions of our methodology to continuous and multinomial cases, as well as to other causal parameters.
In both simulations and an empirical illustration using the National Supported Work Demonstration data, we provide evidence that our procedure performs well compared to existing Bayesian and frequentist approaches. In our Monte Carlo simulations, we find that our method results in improved empirical coverage probabilities, while maintaining very competitive lengths for confidence intervals. This finite sample advantage is also observed over Bayesian methods that rely solely on prior corrections. In particular, we note that our approach leads to more accurate uncertainty quantification and is less sensitive to estimated propensity scores being close to boundary values.
The BvM theorem for parametric Bayesian models is well-established; see, for instance, van1998asymptotic. Its semiparametric version is still being studied very actively when nonparametric priors are used castillo2012gaussian,castillo2015bvm,ray2020causal. To the best of our knowledge, our new semiparametric BvM theorem is the first one that possesses the double robustness property. Our paper is also connected to another active research area concerning Bayesian inference for parameters in econometric models, which is robust to partial or weak identification chen2018MC,giacomini2020robust,andrews2022gmm. The framework and the approach we take is different. Nonetheless, they share the same scope of tailoring the Bayesian inference procedure to new challenges in contemporary econometrics.
This section provides the main setup of the average treatment effect (ATE). We motivate the new Bayesian methodology and detail the practical implementation.
We consider a family of probability distributions $\{P_\eta:\eta\in\mathcal H\}$ for some parameter space $\mathcal H$, where the (possibly infinite dimensional) parameter $\eta$ characterizes the probability model. Let $\eta_0$ be the true value of the parameter and denote $P_0=P_{\eta_0}$, which corresponds to the frequentist distribution generating the observed data.
For individual $i$, consider a treatment indicator $D_i\in\{0,1\}$. The observed outcome $Y_i$ is determined by $Y_i =D_i Y_i(1) +(1 - D_i)Y_i(0)$ where $(Y_i(1), Y_i(0))$ are the potential outcomes of individual $i$ associated with $D_i=1$ or $0$. We now focus on the binary outcome case where both $Y_i(1)$ and $Y_i(0)$ take values in $\{0,1\}$. An extension to multinomial or continuous outcomes is provided in Section (ref). The covariates for individual $i$ are denoted by $X_i$, a vector of dimension $p$, with the distribution $F_0$ and the density $f_0$.\footnote{If $X_i$ does not have a density we can simply consider the conditional density of $(Y_i, D_i)$ given $X_i=x$ instead of the joint density of $(Y_i, D_i, X_i)$.} Let $\pi_0(x)=P_0(D_i=1|X_i=x)$ denote the propensity score and $m_0(d,x)= P_0(Y_i=1|D_i=d, X_i=x)$ the conditional mean. Suppose that the researcher observes independent and identically distributed (i.i.d.) observations of $Z_i=(Y_i,D_i,X_i^\top)^\top$ for $i=1,\dots,n$. The joint density of $Z_i$ is given by $p_{\pi_0,m_0,f_0}$ where
The parameter of interest is the ATE given by $\tau_0=\mathbb{E}_0[Y_i(1)-Y_i(0)]$, where $\mathbb{E}_0[\cdot]$ denotes the expectation under $P_0$. For its identification, we impose the following standard assumption of unconfoundedness and overlap rosenbaum1984,imbens2004,imbens2015causal.
We introduce additional notations from the Bayesian perspective, following the similar setup from ray2020causal. For the purpose of assigning prior distributions to $(\pi, m)$ in the Bayesian procedure, it is convenient to transform them by a link function. We make use of the Logistic function $\Psi(t)=1/(1+e^{-t})$ here. Specifically, we consider the reparametrization of $(\pi, m, f)$ given by $\eta=(\eta^{\pi},\eta^m,\eta^f)$. We index the probability model as $P_{\eta}$, in line with the notation introduced at the first paragraph of this section, where
Below, we write $m_\eta=\Psi(\eta^m)$, $\pi_\eta=\Psi(\eta^\pi)$, and $f_\eta=\exp(\eta^f)$ to make the dependence on $\eta$ explicit. Given any prior on the triplet $(\eta ^{\pi},\eta ^{m},\eta ^{f})$, Bayesian inference on the ATE is achieved by deriving the posterior distribution of
where $\mathbb E_\eta[\cdot]$ denotes the expectation under $P_\eta$. Our aim is to examine large-sample behavior of the posterior of $\tau_{\eta}$ under the true probability distribution $P_0$. In the same vein, the true parameter of interest becomes $\tau_0=\tau_{\eta_0}$.
The construction of our double robust Bayesian procedure in Section (ref) has fundamental connection to the efficient influence function. For any parameter $\eta$, the efficient influence function (hahn1998role,hirano2003efficient) is
for the Riesz representer $\gamma_\eta$, which is given by
We write $\widetilde{\tau}_{0}=\widetilde{\tau}_{\eta_0}$ and $\gamma_{0}=\gamma_{\eta_0}$. Both the prior adjustment and posterior correction of our approach require a pilot estimator for $\gamma_{0}$. Under Assumption (ref), the true Riesz representer $\gamma_0$ is well defined.
We build upon the ATE expression in (ref) to develop our doubly robust inference procedure. Our approach is based on nonparametric prior processes for $\eta^m$ and $\eta^f$. For the latter, we consider the Dirichlet process, which is a default prior on spaces of probability measures. This choice is also convenient for posterior computation via the Bayesian bootstrap; see Remark (ref). For the former, we make use of Gaussian process priors, along with an adjustment that involves a preliminary estimator of $\gamma_0$. Gaussian process priors are also closely related to spline smoothing, as discussed in wahba1990spline. Their posterior contraction properties (see ghosal2017fundamentals), together with excellent finite sample behavior (see rassmusen2006gaussian), make Gaussian process priors popular in the related literature. Since $\tau_{\eta}$ does not depend on $\eta^{\pi}$, the specification of a prior on the propensity score is not required.
We consider pilot estimators $\widehat\pi$ of the propensity score $\pi_0$ and $\widehat m$ of the conditional mean function $m_0$, which both are based on an auxiliary sample. We consider a plug-in estimator for the Riesz representer $\gamma_0$ given by
Below, let $\Gamma_n$ denote the sample average of the absolute value of $\widehat{\gamma}$, which we use for scale normalization in our prior adjustment (see Section (ref) for details). The use of an auxiliary data for pilot estimators simplifies the technical analysis related to the propensity score adjusted priors; see ray2020causal. Also, it provides an effective way to control some negligible higher-order terms, see our Lemma (ref) in the Supplemental Material and the related discussion about the sample splitting in the DML type methods on Page C6 of chernozhukov2018double. In practice, we use the full data twice and do not split the sample, as we have not observed any over-fitting or loss of coverage thereby. Algorithm (ref) describes our double robust Bayesian inference procedure.
Given the draws from the corrected posterior calculated in Algorithm (ref), we obtain the point estimate and credible set as follows. The Bayesian point estimator is $\overline{\tau}_{\eta}=\frac{1}{S}\sum_{s=1}^S \check{\tau}_\eta^{s}$. The $100\cdot(1-\alpha)\%$ credible set for the ATE parameter $\tau_0$ is given by
where $q_n(a)$ denotes the $a$-th quantile of $\{\check{\tau}_\eta^{s}:s=1,\ldots,S\}$.
For the implementation of our pilot estimator $\widehat\gamma$ given in (ref), we recommend using propensity scores estimated by the Logistic Lasso. For the implementation of the pilot estimator $\widehat m$, we adopt the posterior mean of $m_{\eta}$ generated from a Gaussian process prior without adjustment, as in ghosal2006binary. Section (ref) provides more implementation details. To approximate the posterior distribution, we make use of the Laplace approximation, but one can also resort to the Markov Chain Monte Carlo (MCMC) algorithms. The parameter $\sigma_n$ controls the relative weight placed on the prior adjustment relative to the standard unadjusted prior on $\eta^m$ (e.g., a Gaussian prior with a squared exponential covariance function). Regarding the tuning parameter $\sigma_n$, we emphasize that our finite sample results are not sensitive to its choice, as we show in Supplemental Appendix (ref).
In this section, we derive the Bernstein-von Mises (BvM) theorem which establishes the asymptotic equivalance between our Bayesian procedure and the frequentist-type semiparametric efficient one for the ATE. We consider an asymptotically efficient estimator $\widehat{\tau}$ with the following linear representation:
where $\widetilde{\tau}_0=\widetilde{\tau}_{\eta_0}$ is the efficient influence function in accordance with (ref). Below, we denote $Z^{(n)}=(Z_1,\ldots, Z_n)$. By virtue of the BvM Theorem, two conditional distributions $\sqrt{n}(\tau_{\eta}-\widehat{\tau})|Z^{(n)}$ and $\sqrt{n}(\widehat{\tau}-\tau_{\eta})|\eta=\eta_0$ are asymptotically equivalent under the underlying sampling distribution. Another important consequence of the BvM theorem is about the asymptotic normality and efficiency of the Bayesian point estimator. That is, $\sqrt n(\overline\tau_{\eta}-\tau_0)$ is asymptotically normal with mean zero and variance $\textsc v_0=\mathbb{E}_0\left[ \widetilde{\tau}_0^2(Z_i)\right]$. Thus, $\overline\tau_{\eta}$ achieves the semiparametric efficiency bound of hahn1998role.
Our prior correction through the Riesz representer $\gamma_0$ is motivated by the least favorable direction of Bayesian submodels. We first provide such least favorable calculations, which are closely linked to the semiparametric efficiency. Consider the one-dimensional submodel $t\mapsto \eta_t$ defined by the path
for a given direction $(\mathfrak{p}, \mathfrak{m},\mathfrak{f})$ with $\int \mathfrak{f}(x)f(x)\,\mathrm{d}x=0$. The difficulty of estimating the parameter $\tau_{\eta_{t}}$ for the submodels depends on the direction $(\mathfrak{p}, \mathfrak{m},\mathfrak{f})$. Among them, let $ \xi_{\eta}= (\xi_{\eta}^{\pi},\xi_{\eta}^{m},\xi_{\eta}^{f})$ be the least favorable direction that is associated with the most difficult submodel. It yields the largest asymptotic optimal variance for estimating $\tau_{\eta_{t}}$ among all submodels. Let $p_{\eta_t}$ denote the joint density of $Z$ depending on $\eta_t:= (\pi_t, m_t, f_t)$. Taking derivative of the logarithmic density $\log p_{\eta_t} (z)$ with respect to $t$ and evaluating at $t=0$ gives the score operator:
where $B_{\eta}^{\pi}\mathfrak{p}(z) = (d-\pi_\eta(x))\mathfrak{p}(x)$, $B_{\eta}^{m}\mathfrak{m}(z)= (y-m_\eta(d,x))\mathfrak{m}(d,x)$ and $B_{\eta}^{f}\mathfrak{f}(z)= \mathfrak{f}(x)$. The least favorable direction is defined as the solution $\xi_\eta$ which solves the equation $B_{\eta}\xi_{\eta}=\widetilde{\tau}_{\eta}$, see ghosal2017fundamentals. We immediately obtain the following.
Lemma (ref) motivates the adjustment of the prior distribution as considered in our Bayesian procedure in Section (ref). Our prior correction, which takes the form of the (estimated) least favorable direction, provides an exact invariance under a shift of nonparametric components in this direction. It provides additional robustness against posterior inaccuracy in the “most difficult direction”, i.e., the one inducing the largest bias in the average treatment effects. We also note that Lemma (ref) extends the result in Section 2.1 in ray2020causal for the missing data problem, which is equivalent as observing only one arm (either the treatment or control arm), to the context of ATE estimation that involves both arms.
We now provide additional notations and assumptions. The posterior distribution plays an important role in the following analysis and is given by
where $p_{\pi,m}$ denotes the conditional density of $(Y_i, D_i) $ given $X_i$, given by (ref) divided by the marginal density of $X_i$. We write $\mathcal{L}_{\Pi}(\sqrt{n}(\tau_\eta-\widehat{\tau})|Z^{(n)})$ for the marginal posterior distribution of $\sqrt{n}(\tau_\eta-\widehat{\tau})$. We focus on the case that $\eta^{\pi }$ has a prior that is independent of the prior for $(\eta^m,F)$. Because the likelihood function ((ref)) factorizes into $(\eta^m,\eta^{\pi},F)$ separately, the posterior of $\eta^{\pi }$ is also independent of the posterior for $(\eta^m,F)$. Due to the fact that $\tau_{\eta}$ does not depend on $\eta^{\pi}$, it is unnecessary to further discuss a prior or posterior distribution on $\eta^{\pi}$.
We first introduce high-level assumptions and discuss primitive conditions for those in the next section. Below, we consider some measurable sets $\mathcal H^m_n$ of functions $\eta^m$ such that $\Pi(\eta^m\in\mathcal{H}^m_n|Z^{(n)})\to_{P_0} 1$. W e also denote $\mathcal{H}_n=\{\eta:\eta^m\in\mathcal{H}_n^m\}$ when we index the conditional mean function $m_{\eta}$ by its subscript $\eta$. We introduce the notation $\|\phi\|_{2, F_0}:= \sqrt{\int \phi^2(x)\,\mathrm{d}F_0(x)}$ for all $\phi\in L^2(F_0):=\{\phi:\|\phi\|_{2, F_0}<\infty\}$, as well as the supremum norm $\|\cdot\|_\infty$. For two sequences $\{a_n\}$ and $\{b_n\}$ of positive numbers, we write $a_n \lesssim b_n$ if $\limsup_{n\to\infty} (a_n / b_n)<\infty$, and $a_n \sim b_n$ if $a_n \lesssim b_n$ and $b_n \lesssim a_n$.
We adopt the standard empirical process notations as follows. For a function $h$ of a random vector $Z_i$ that follows distribution $P_0$, we let $P_0[h]=\int h(z)\,\mathrm{d}P(z),\mathbb{P}_n[h]=n^{-1}\sum_{i=1}^{n}h(Z_i)$, and $\mathbb{G}_n[h]=\sqrt n\left(\mathbb{P}_n-P_0\right)[h]$. Below, we make use of the notations $\bar{m}_{\eta}(\cdot)=m_{\eta}(1,\cdot)-m_{\eta}(0,\cdot)$ and $\bar{m}_{0}(\cdot)=m_{0}(1,\cdot)-m_{0}(0,\cdot)$.
Recall the propensity score-dependent prior on $m$ given by $m_\eta(\cdot) = \Psi\left(\eta^m(\cdot)\right)$ where $\eta^m(\cdot )=W^m(\cdot) + \lambda\,\widehat \gamma(\cdot)$. The restriction on $\lambda$ is made through its hyperparameter $\sigma_n>0$.
Discussion of Assumptions: Assumption (ref) imposes sufficiently fast convergence rates for the pilot estimators for the conditional mean function $m_{0}$ and the propensity score $\pi _{0}$. When considering frequentist pilot estimators, these rate conditions can be justified by adopting the recent proposals of ChernozhukovNeweySingh2020a, ChernozhukovNeweySingh2020b. One can also use Bayesian point estimators such as the posterior mean of the Gaussian process for $\widehat{m}$ and $\widehat{\pi}$. The posterior convergence rate for the conditional mean $m_{\eta}$ can be derived in the same spirit of ray2020causal. The rate conditions in Assumption 2 also resemble conditions (i) and (ii) of Theorem 1 of farrell2015 in the context of frequentist estimation. Remark (ref) illustrates that under classical smoothness assumptions, this assumption is less restrictive than the method of ray2020causal or other approaches for semiparametric estimation of ATEs as found in chen2008semiparametric or farrell2021deep. Assumption (ref) incorporates Conditions (3.9) and (3.10) from Theorem 2 in ray2020causal, and it is imposed to check the invariance property of the adjusted prior distribution. These restrictions are mild and extend beyond the Gaussian processes considered in Section (ref) for concreteness.
Assumption (ref) restricts the functional class $\mathcal{G}_{n}$ to form a $P_{0}$-Glivenko--Cantelli class; see Section 2.4 of van1996empirical. This imposes a new stochastic equicontinuity condition, as ((ref)) restricts a product structure involving $\widehat\gamma$ and $m_{\eta}$, which further relaxes the corresponding condition from ray2020causal, namely, $\sup_{\eta \in \mathcal{H}^{m}_{n}}\mathbb{G}_{n} [m_{\eta}-m_{0}] = o_{P_{0}}(1)$. In the next section, we demonstrate that our formulation allows for double robustness under H\"older classes (see Remark (ref)). Hence, the complexity of the functional class $(m_{\eta}-m_{0})$ can be compensated by sufficient regularity of the corresponding Riesz representer and vice versa. A condition similar to our Assumption (ref) is also used in the frequentist literature; see Section 2 of Benkeser, Carone, van der Laan, and Gilbert (benkeser2017doubly). Nonetheless, the technical argument differs substantially from the frequentist's study, because we mainly need the condition ((ref)) to control changes in the likelihood under perturbations along the estimated and true least favorable directions. This is unique to Bayesian analysis with nonparametric priors.
We now establish a new Bernstein–von Mises theorem, which establishes the asymptotic normality of the posterior distribution, modulo a “bias term”. In a next step, we show that posterior correction, as proposed in our procedure, eliminates this “bias term”. This asymptotic equivalence result is established using the bounded Lipschitz distance. For two probability measures $P,Q$ defined on a metric space $\mathcal{Z}$, we define the bounded Lipschitz distance as
where
Here, $\|\cdot\|_{\ell_2}$ denotes the vector $\ell_2$ norm.
Below is our main statement about the asymptotic behavior of the posterior distribution of $\tau_{\eta}$. As in the modern Bayesian paradigm, the exact posterior is rarely of closed-form, and one needs to rely on certain Monte Carlo simulations, such as the implementation procedure in Section (ref), to approximate this posterior distribution, as well as the resulting point estimator and credible set.
We emphasize that the above BvM theorem is not feasible for applications, because it depends on the “bias term” $b_{0,\eta}$, which depends on the unknown conditional mean $m_0$. Nonetheless, it provides an important theoretical benchmark. One can follow the existing literature on semiparametric BvM theorems to impose the so-called “no-bias" condition, but this generally leads to strong smoothness restrictions and may not be satisfied when the dimensionality of covariates is large relative to the smoothness properties of the underlying functions; see the discussion on page 395 of van1998asymptotic.
This “bias term” in our context consists of two key components, with the first involving unknown true functions and the second depending on the posterior of $m_{\eta}$. We consider pilot estimators for the unknown functional parameters in $b_{0,\eta}$. The correction term $\widehat{b}_{\eta}$, as introduced in (ref), results in a feasible Bayesian procedure that satisfies the BvM theorem under double robustness, as demonstrated below.
We now show how Theorem (ref) can provide frequentist justification of Bayesian methods to construct the point estimator and the confidence sets. Recall that $\overline{\tau}_{\eta}$ represents the posterior mean. Introduce a Bayesian credible set $\mathcal{C}_n(\alpha)$ for $\tau_\eta$, which satisfies $\Pi(\tau_\eta\in \mathcal{C}_n(\alpha)|Z^{(n)})=1-\alpha$ for a given nominal level $\alpha\in(0,1)$. The next result shows that $\mathcal{C}_n(\alpha)$ also forms a confidence interval in the frequentist sense for the ATE parameter whose coverage probability under $P_0$ converges to $1-\alpha$.
To the best of our knowledge, this is the first BvM theorem that entails the double robustness. We discuss the distinction from Theorem 2 in ray2020causal. Their work laid the theoretical foundation for Bayesian inference based on propensity score adjusted priors. Specifically, under this prior adjustment, they establish a BvM result under weak regularity conditions on the propensity score function. Our analysis differs from ray2020causal in two crucial ways. First, we improve on their Lemma 3 by showing that it is possible to verify the prior stability condition for propensity score-adjusted priors under the product structure in Assumption (ref), modulo the “bias term" $b_{0,\eta}$. This separation is essential to identify the source of the restrictive condition, such as the Donsker property on $m_{\eta}$, which is mainly used to eliminate $b_{0,\eta}$. Second, our proposal introduces an explicit debiasing step, borrowing key insights from recent developments in the DML literature.
We illustrate the general methodology by placing a particular Gaussian process prior on $\eta^m(d,\cdot)$ in relation to the conditional mean functions for $d\in\{0,1\}$. The Gaussian process regression has been extensively used among the machine learning community, and started to gain popularity among economists kasy2018tax. We provide primitive conditions used in our main results in the previous section. In addition, we provide details on the implementation using Gaussian process priors and discuss the data-driven choices of tuning parameters.
Let $(W(t):t\in\mathbb{R}^p)$ be a generic centered and homogeneous Gaussian random field with covariance function of the following form $\mathbb{E}[W(s)W(t)]=\phi(s-t)$, for a given continuous function $\phi:\mathbb{R}^p\mapsto \mathbb{R}$. We consider $W(t)$ as a Borel measurable map in the space of continuous functions on $[0,1]^p$, equipped with the supremum norm $\Vert \cdot\Vert_{\infty}$. The Gaussian process is completely determined by the covariance function. For example, the covariance function of the squared exponential process is given by $\mathbb{E}[W(s)W(t)]=\exp(-\Vert s-t\Vert_{\ell_2}^2)$, as its name suggests. In this section, we focus on the squared exponential process prior, which is one of the most commonly used priors in applications; see rassmusen2006gaussian and murphy2023pml. We also consider a rescaled Gaussian process $\big(W(a_nt):\,t\in [0,1]^p\big)$ for some $a_n>0$. Intuitively speaking, $a_n^{-1}$ can be thought as a bandwidth parameter. For a large $a_n$ (or equivalently a small bandwidth), the prior sample path $t\mapsto W(a_nt)$ is obtained by shrinking the long sample path $t\mapsto W(t)$. Thus, it employs more randomness and becomes suitable as a prior model for less regular functions, see van2008gaussian,van2009adaptive.
Below, $\mathcal{C}^{s_m}([0,1]^p)$ denotes a H\"older space with the smoothness index $s_m$. Specifically, we illustrate our theory with the case where $m_0(d,\cdot)\in \mathcal{C}^{s_m}([0,1]^{p})$ for $d\in\{0,1\}$. Given such a H\"older-type smoothness condition, we choose
Under ((ref)), a rescaled Gaussian process $\big(W(a_nt):\,t\in [0,1]^p\big)$ induces the posterior contraction rate for the conditional mean function $m_\eta(d,\cdot)$ to be $\varepsilon_n=n^{-s_m/(2s_m+p)}(\log n)^{s_m(1+p)/(2s_m+p)}$; see Section 11.5 of ghosal2017fundamentals. The particular choice of $a_n$ mimics the corresponding kernel bandwidth based on kernel smoothing methods. Other choices of $a_n$ will generally make the convergence rate slower. Nonetheless, as long as the propensity score is estimated with a sufficiently fast rate, our BvM theorem still holds. The next proposition illustrates our general theory when we adopt the rescaled squared exponential process prior for the conditional mean function. We use the superscript $m$ for the prior process $W^m$ to signify this relationship.
We provide details on the Gaussian process prior placed on $\eta^m(d,x)$ and its posterior computation. Algorithm (ref) sets the adjusted prior as $\eta^m(d,x)=W^m(d,x) + \lambda\,\widehat \gamma(d,x)$: In our implementation, we choose the first component $W^m(d,x)$ to be a zero-mean Gaussian process with the commonly used squared exponential covariance function rassmusen2006gaussian. That is, $K\left((d,x),(d',x')\right):= \nu^2 \exp\left(-a_{0n}^{2}(d-d^\prime)^2/2-\sum_{l=1}^{p}a_{ln}^{2}(x_{l}-x^\prime_{l})^2/2\right)$ where the hyperparameter $\nu^2$ is the kernel variance and $a_{0n},\ldots,a_{pn}$ are rescaling parameters that reflect the relevance of treatment and each covariate in predicting $\eta^m$. They are selected by maximizing the marginal likelihood. Conditional on the data used to obtain the propensity score estimator $\widehat\pi$, the prior for $\eta^m$ has zero mean and the covariance kernel $K^c$, which includes an additional term based on the estimated Riesz representer $\widehat\gamma$. It is given by $K^c\left((d,x),(d^\prime,x^\prime)\right) = K\left((d,x),(d^\prime,x^\prime)\right) + \sigma_n^2\widehat \gamma(d,x)\,\widehat \gamma(d^\prime,x^\prime),$ cf. related constructions from ray2019debiased and ray2020causal. The parameter $\sigma_n$, representing the standard deviation of $\lambda$, controls the weight of the prior adjustment relative to the standard Gaussian process. The choice $\sigma_n=(\log n)/(\sqrt{n}\,\Gamma_n)$, where $\Gamma_n= n^{-1}\sum_{i=1}^n\vert\widehat\gamma(D_i,X_i)\vert$, as specified in Algorithm (ref), satisfies the conditions $\sigma_n\lesssim 1$ and $n\sigma^2_{n}\to\infty$ in Assumption (ref), with probability approaching one. It is similar to the choice suggested by ray2019debiased, where $\sigma_n$ is proportional to $1/(\sqrt{n}\,\Gamma_n)$. The factor $\Gamma_n$ normalizes the second term (adjustment term) of $K^c$ to have the same scale as the unadjusted covariance $K$. Supplemental Appendix (ref) shows that the finite sample performance of the double robust Bayesian approaches remains stable across different choices of $\sigma_n$.
Utilizing Gaussian process priors with zero mean and covariance function $K^c$, and incorporating the available data, we generate posterior draws of the vector $ \left[\eta^m(d,X_1),\cdots,\eta^m(d,X_n)\right]^{\top}$ for $d\in \{0,1\}$. This can be achieved through the Laplace approximation method detailed in Supplemental Appendix (ref).
For the implementation of the pilot estimator $\widehat\gamma$ given in (ref), we recommend Logistic Lasso for the propensity score, with the penalty parameter chosen by cross-validation friedman2010regularization. As a pilot estimator $\widehat{m}$ in Algorithm (ref) for posterior correction, we use the uncorrected posterior mean $\sum_{s=1}^S m_\eta^{s}/S$, where $m_\eta^{s}$ is calculated following Step (a) of posterior computation in Algorithm (ref), but with a Gaussian process prior without adjustment, i.e., $\Psi\left(W^m(d,\cdot)\right)$. When the rescaling parameter $a_n$ is as stated in Proposition (ref), the convergence rate of $\widehat{m}$ is $O_{P_0}\big((n/\log n)^{-s_m/(2s_m+p)}\big)$. This can be shown by combining Theorems 11.22, 11.55 and 8.8 from ghosal2017fundamentals.
In this section, we apply our method to one version of the Lalonde--Dehejia--Wahba data that contains a treated sample of 185 men from the National Supported Work (NSW) experiment and a control sample of 2490 men from the Panel Study of Income Dynamics (PSID). The data has been used by lalonde1986, dehejia1999causal, abadie2011bias, and armstrong2021finite, among others. We refer readers to lalonde1986, and dehejia1999causal for reviews of the data.\footnote{The data is available on Dehejia's website: \href{http://users.nber.org/ rdehejia/nswdata2.html}{http://users.nber.org//$\sim$ rdehejia/nswdata2.html}.}
In this section, we consider a simulation study where the observations are randomly drawn from a large sample generated by the Wasserstein Generative Adversarial Networks (WGAN) method from the the Lalonde--Dehejia--Wahba data, see athey2021using. We view their simulated data as the population and repeatedly draw our simulation samples (each consisting of 185 treated observations and 2490 control observations) for each of the $1000$ Monte Carlo replications. We slightly depart from previous studies by focusing on a binary outcome $Y$: the employment indicator for the year 1978, which is defined as an indicator for positive earnings. The treatment $D$ is the participation in the NSW program. We are interested in the average treatment effect of the NSW program on the employment status. For the set of covariates, we follow abadie2011bias and include nine variables: age, education, black, Hispanic, married, earnings in 1974, earnings in 1975, unemployed in 1974, and unemployed in 1975. We implement our double robust Bayesian method (DR Bayes) following Algorithm (ref), using $S=5000$ posterior draws and the pilot estimator $\widehat \gamma$ and $\widehat m$, as detailed at the end of Section (ref). We compare DR Bayes to two other Bayesian procedures: First, we consider the prior adjusted Bayesian method (PA Bayes) proposed by ray2020causal, which constructs the point estimate and credible interval based on $\tau_{\eta}^s$ in ((ref)). Second, we examine an unadjusted Bayesian method (Bayes) which is also based on $\tau_{\eta}^s$ but is generated using Gaussian process priors without the adjustment.
We also compare our method to frequentist estimators. Match/Match BC corresponds to the nearest neighbor matching estimator and its bias-corrected version by abadie2011bias, which adjusts for differences in covariate values through regression. DR TMLE corresponds to the doubly robust targeted maximum likelihood estimator by benkeser2017doubly. DML refers to the double/debiased machine learning estimator from chernozhukov2017double, where the nuisance functions $\pi_0$ and $m_0$ are estimated using random forests (which outperformed DML combined with other nuisance function estimators, such as Lasso, in our simulation setup). Since the job-training data contains a sizable proportion of units with propensity score estimates very close to $0$ and $1$, we follow crump2009dealing and discard observations with the estimated propensity score outside the range $[t, 1-t]$, with the trimming threshold $t\in\{0.10, 0.05, 0.01\}$.\footnote{crump2009dealing suggested a simple rule of thumb with a threshold of $t=0.10$, while athey2021using used $t=0.05$. Applying the optimal trimming rule proposed by crump2009dealing to our simulated samples yields an average optimal trimming threshold $0.073$.}
Table (ref) presents the finite sample performance of the Bayesian and frequentist methods mentioned above. We use the full data twice in computing the prior/posterior adjustments and the posterior distribution of the conditional mean function. Supplemental Appendix (ref) reports the performance of DR Bayes using sample splitting, which results in similar coverage but a larger credible interval length due to the halved sample size.
Concerning the Bayesian methods for estimating the ATE, Table (ref) reveals that unadjusted Bayes yields highly inaccurate coverage except for the case with trimming constant $t=0.01$. If the prior is corrected using the propensity score adjustment, the results improve significantly. Nevertheless, our DR Bayes method demonstrates two further improvements: First, DR Bayes leads to smaller average confidence lengths in each case while simultaneously improving the coverage probability. This can be attributed to a reduction in bias and/or more accurate uncertainty quantification via our posterior correction. Second, when the trimming threshold is small (i.e., $t=0.01$), propensity score estimators can be less accurate, leading to reduced coverage probabilities of PA Bayes. Our double robust Bayesian method, on the other hand, is still able to provide accurate coverage probabilities. In other words, DR Bayes exhibits more stable performance than PA Bayes with respect to the trimming threshold.\footnote{In additional simulations without trimming ($t=0$), we find that all double robust methods, including DR Bayes, substantially under-cover and/or inflate the length of their confidence intervals. This is consistent with crump2009dealing, who point out that propensity score estimates close to the boundaries tend to induce substantial bias and large variances in estimating the ATE. We also note that unadjusted Bayes severely undercovers in this case.}
Our DR Bayes also exhibits encouraging performances when compared to frequentist methods. It provides a more accurate coverage than bias-corrected matching, DR TMLE and DML. Compared with the matching estimator that exhibits a similarly good coverage performance, DR Bayes yields considerably shorter credible intervals.
We apply the Bayesian and frequentist methods considered above to the Lalonde--Dehejia--Wahba data. Similar to the simulation exercise, we consider a varying choice of the threshold $t\in \{0.10,0.05,0.01\}$.\footnote{Applying the optimal trimming rule proposed by crump2009dealing yields an optimal threshold of $0.064$.} The ATE point estimates and confidence intervals are presented in Table (ref). As a benchmark, the experimental data that uses both treated and control groups in NSW ($n=445$) yields an ATE estimate (treated-control mean difference) of $0.111$ with a $95\%$ confidence interval $[0.026, 0.196]$.
As we see from Table (ref), the unadjusted Bayesian method yields larger estimates. The adjusted Bayesian methods (PA and DR Bayes), on the other hand, produce estimates comparable to the experimental estimate. PA Bayes finds that the job training program enhanced the employment by $9.0\%$ to $17.0\%$ across different trimming thresholds, and DR Bayes estimates the effect from $12.1\%$ to $18.4\%$. Among frequentist estimators, the matching estimator and its bias-corrected version produce similar estimates as PA and DR Bayes, but with wider confidence intervals. DR TMLE produces negative estimates for $t=0.10$ when all other estimates are positive. For $t=0.10$ and $0.05$, DML yields similar point estimates as PA and DR Bayes, but with less estimation precision. In the case $t=0.01$ where the overlapping condition is closer to violation, however, its point estimate and confidence interval length become considerably larger than other methods.
This section extends the binary variable $Y$ to encompass general cases, including continuous, counting, and multinomial outcomes. First, we examine the class of single-parameter exponential families, where the conditional density function is solely determined by the nonparmatric conditional mean function. This covers continuous outcomes and counting variables. Second, we consider the “vector" case of exponential families for multinomial outcomes. For both classes, we derive the novel correction to the Bayesian procedure and delegate more technical discussions to Supplemental Appendices (ref) and (ref). Additionally, we outline extensions to other causal parameters of interest.
In this part, we assume that the distribution of $Y_i$ conditional on $D_i$ and $X_i$ belongs to the “single-parameter" exponential family, where the unknown parameter is the nonparametric conditional mean function $m(d,x)=\mathbb{E}[Y_i|D_i=d,X_i=x]$. The conditional density function is given by
where $A(m)= \log\int c(y)\exp\left[q(m) ay\right]\mathrm{d}y$, and the function $q(\cdot)$ links the mean to the “natural parameter” of the exponential family. We also restrict the sufficient statistic to be linear in $y$.
The family ((ref)) not only encompasses the Bernoulli distribution (with $q(m)=\log(m/(1-m))$, $A(m)=-\log(1-m)$, and $c(y)=a=1$), as considered in the previous sections, but also allows for counting and continuous outcomes. For instance, when $a=1$, the Poisson distribution corresponds to the choices $c(y)= 1/(y!)$, $q(m)=\log m$, and $A(m)=m$, while the exponential distribution is represented by $c(y)=1$, $q(m)=-1/m$, and $A(m)=\log m$. Furthermore, the normal distribution with $\text{Var}(Y|D,X)=\sigma^2$ for some $\sigma>0$, is captured by $c(y)=\exp(-y^2/(2\sigma^2))/\sqrt{2\pi\sigma^2}$, $q(m)=m/\sigma$, $A(m)=m^2 / (2\sigma^2)$, and $a=1/\sigma$. We emphasize that model (ref) does not impose functional form assumptions on the conditional mean function $m$. The joint density of $(Y_i,D_i,X_i)$ can be written as
We consider the same reparametrization of $(\pi, m, f)$ as in (ref) except that now the second component of $\eta$ uses the general link function $q$ satisfying $\eta^{m} = q(m)$. We now state the least favorable direction for the exponential family case, which serves as motivation for the prior adjustment.
For the outcome family with $a=1$, which includes Bernoulli, Poisson and exponential distributions, the least favorable direction for ATE estimation coincides with the one as given in Lemma (ref). To implement the double robust Bayesian procedure for general outcomes, one can still follow Algorithm (ref), with the logistic function $\Psi$ replaced by the inverse link function $q^{-1}$. For the normal (homoscedastic) outcome where prior adjustment $\lambda\widehat\gamma(d,x)$ in Algorithm (ref) becomes $\lambda\widehat\gamma(d,x)/a$, the hyperparameter $a$ can be determined together with other parameters of the Gaussian process by optimizing the marginal likelihood as in ray2019debiased. Proposition (ref) in the Supplemental Material provides primitive conditions for the BvM Theorem to hold under double robust smoothness conditions.
We now assume that the dependent variable $Y_i$ takes values in a finite set, specifically $Y_i\in\{0,1, \dots,J\}$. The ATE can then be written as $\tau_\eta=\sum_{j=0}^J j\, \mathbb E_\eta\left[m_{\eta,j}(1,X) - m_{\eta, j}(0,X) \right]$, where the choice probabilities are $m_{\eta, j}(d,x) = \Psi_j\left(\eta^{m_1},\cdots,\eta^{m_J}\right)(d,x)$ with the multinomial logit specification:
for $j=1,\ldots, J.$ The multinomial logit specification implies $m_{\eta, 0}(d,x) =1-\sum_{j=1}^{J}m_{\eta, j}(d,x)$. We now provide the least favorable direction for multinomial outcomes in the presence of multinomial outcomes and discuss its consequences for prior adjustment below.
We emphasize that the least favorable direction calculation is not a trivial extension of hahn1998role or ray2020causal. This is because there are $J$ nonparametric components involved in the conditional probability function of the multinomial outcomes given covariates, and we need to consider the perturbation of those $J$ components together. Nonetheless, we show that the efficient influence function is of the same generic form as derived in hahn1998role. In the proof of Lemma (ref), we compute the derivative of the parameter mapping along the path considered herein. We derive inner products involving the least favorable direction for each nonparametric component consisting of the conditional choice probabilities. The extension to the multinomial case had not been considered in the literature to our knowledge, and it offers a result of independent interest.
Lemma (ref) motivates the following modification of our double robust Bayesian estimator based on the propensity score-dependent prior on $m_{\eta,j}$ for $1\leq j\leq J$:
where $W^{m_j}(d,\cdot)$ is a continuous stochastic process independent $\lambda\sim N(0,\sigma_n^2)$ for $\sigma_n>0$. We may then follow the implementation as described in Section (ref) using $m_\eta(d,x)=\sum_{j=0}^J j\, m_{\eta,j}(d,x)$.
We now extend our procedure to general linear functionals of the conditional mean function. We do so only for binary outcomes, as the modification to other types of outcomes follows as above. Recall that the observable data consists of $i.i.d.$ observations of $Z=(Y,D,X^\top)^\top$. The causal parameter of interest is $\tau_0=\mathbb{E}_0[\psi(Z,m_0)]$, where the function $\psi$ is linear with respect to the conditional mean function $m_0$. We introduce the Riesz representer $\gamma_0(d,x)$ satisfying $\mathbb{E}_0[\psi(Z,m)]=\mathbb{E}_0[\gamma_{0}(D,X)m(D,X)]$. Let $\widehat{m}$ and $\widehat{\gamma}$ be pilot estimators for the conditional mean and Riesz representer, respectively, computed over an auxiliary sample. Our double robust Bayesian procedure can be extended by considering the corrected posterior distribution for $\tau_\eta$ as follows: $ \check{\tau}_\eta^{s}= \sum_{i=1}^n M_{ni}^s \psi(Z_i,m^s_\eta)-n^{-1}\sum_{i=1}^n \boldsymbol{\tau}[m_\eta^s-\widehat m](Z_i)$, $s=1,\ldots,S$, where here $\boldsymbol{\tau}[m](z):=\psi(z,m)+\widehat{\gamma}(d,x)(y-m(d,x))$. The derivations of the least favorable directions in the following two examples are provided in Supplemental Appendix (ref).