The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
68,712 characters
Double Robust Bayesian Inference on Average Treatment Effects
\maketitle
\begin{abstract}
\vskip -.1cm
{\small
We propose a double robust Bayesian inference procedure on the average treatment effect (ATE) under unconfoundedness.
For our new Bayesian approach, we first adjust the prior distributions of the conditional mean functions, and then correct the posterior distribution of the resulting ATE.
Both adjustments make use of pilot estimators motivated by the semiparametric influence function for ATE estimation.
We prove asymptotic equivalence of our Bayesian procedure and efficient frequentist ATE estimators by establishing a new semiparametric Bernstein-von Mises theorem under double robustness; i.e., the lack of smoothness of conditional mean functions can be compensated by high regularity of the propensity score and vice versa.
Consequently, the resulting Bayesian credible sets form confidence intervals with asymptotically exact coverage probability.
In simulations, our method provides precise point estimates of the ATE through the posterior mean and credible intervals that closely align with the nominal coverage probability. Furthermore, our approach achieves a shorter interval length in comparison to existing methods. We illustrate our method in an application to the National Supported Work Demonstration following \cite{lalonde1986} and \cite{dehejia1999causal}.
}
\end{abstract}
\noindent{\footnotesize \noindent \textsc{Keywords:} Average treatment effects, unconfoundedness, double robustness, nonparametric Bayesian inference, Bernstein–von Mises theorem, Gaussian processes. }
\section{Introduction}
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 \cite{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 \cite{robins1997} or \cite{ritov2014}.
Once the prior distribution was adjusted using some pre-estimated propensity score, \cite{ray2020causal} recently established a novel semiparametric BvM theorem under weaker smoothness requirement for the propensity score function.\footnote{Strictly speaking, the main objective in \cite{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 \cite{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 \cite{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 \cite{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 \cite{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 \cite{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, \cite{van1998asymptotic}. Its semiparametric version is still being studied very actively when nonparametric priors are used \citep{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 \citep{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.
\section{Setup and Implementation}\label{sec:model}
This section provides the main setup of the average treatment effect (ATE). We motivate the new Bayesian methodology and detail the practical implementation.
\subsection{Setup}
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{sec:extension}.
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
\begin{equation}\label{x_den}
p_{\pi,m,f}(z)=\pi(x)^d (1-\pi(x))^{1-d}m(d,x)^{y} (1-m(d,x))^{(1-y)}f(x).
\end{equation}
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 \citep{rosenbaum1984,imbens2004,imbens2015causal}.
\begin{assumption}\label{Ass:unconfounded}
(i) $(Y_i(0),Y_i(1)) ~ \raisebox{0.05em}{\rotatebox[origin=c]{90}{$\models$}} ~ D_i \mid X_i$ and (ii) there exists $\bar\pi>0$ such that $\bar\pi <\pi_0(x)< 1-\bar\pi$ for all $x$ in the support of $F_0$.
\end{assumption}
We introduce additional notations from the Bayesian perspective, following the similar setup from \cite{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
\begin{eqnarray}\label{repar}
\eta^{\pi} = \Psi^{-1}(\pi), ~~\eta^{m} = \Psi^{-1}(m),~~\eta^f =\log f.
\end{eqnarray}
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
\begin{equation}\label{ate}
\tau_\eta=\mathbb E_\eta\left[m_\eta(1,X)-m_\eta(0,X)\right],
\end{equation}
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{sec:method_outline} has fundamental connection to the efficient influence function. For any parameter $\eta$, the efficient influence function (\cite{hahn1998role,hirano2003efficient}) is
\begin{align}\label{eif_ate}
\widetilde{\tau}_{\eta}(z) &= m_\eta(1,x)-m_\eta(0,x)+ \gamma_\eta(d,x)(y-m_\eta(d,x))-\tau_\eta
\end{align}
for the Riesz representer $\gamma_\eta$, which is given by
\begin{equation}\label{riesz:def}
\gamma_\eta(d,x)=\frac{d}{\pi_\eta(x)}-\frac{1-d}{1-\pi_\eta(x)}.
\end{equation}
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{Ass:unconfounded}, the true Riesz representer $\gamma_0$ is well defined.
\subsection{Double Robust Bayesian Point Estimators and Credible Sets}\label{sec:method_outline}
We build upon the ATE expression in \eqref{ate} 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{rem:BB}. 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 \cite{wahba1990spline}. Their posterior contraction properties (see \cite{ghosal2017fundamentals}), together with excellent finite sample behavior (see \cite{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
\begin{align}\label{Riesz_est}
\widehat{\gamma}(d,x)=\frac{d}{\widehat\pi(x)}-\frac{1-d}{1-\widehat\pi(x)}.
\end{align}
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{sec:implement:GP} for details).
The use of an auxiliary data for pilot estimators simplifies the technical analysis related to the propensity score adjusted priors; see \cite{ray2020causal}. Also, it provides an effective way to control some negligible higher-order terms, see our Lemma \ref{lemma:negligible} in the Supplemental Material and the related discussion about the sample splitting in the DML type methods on Page C6 of \cite{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{algorithm} describes our double robust Bayesian inference procedure.
\begin{algorithm}[H]
\caption{Double Robust Bayesian Procedure}
\label{algorithm}
\begin{algorithmic}
\STATE \textbf{Input:} Data $Z_i=(Y_i,D_i,X_i^\top)^\top$ for $i=1,\dots,n$, number of posterior draws $S$, initial estimators $\widehat \gamma$ and $\widehat m$,
and $\lambda \sim N(0,\sigma_n^2)$ where $\sigma_n=\left(\log n\right)/(\sqrt n\, \Gamma_n)$.
\STATE \textbf{Prior Specification:}
\STATE (a) Select a Gaussian process prior $W^m$.
\STATE
(b) Set an adjusted prior for $m_\eta(d,X_i) = \Psi\left(\eta^m(d,X_i)\right)$, where $ \eta^m(d,X_i)=W^m(d,X_i) + \lambda\,\widehat \gamma(d,X_i)$. \\
\textbf{Posterior Computation:}
\FOR{$s=1,\ldots, S$}
\STATE (a)
Generate the $s$-th draw of the posterior of $(m_\eta(d, X_i))_{i=1}^n$ using the adjusted prior and the data; denote it as $(m^s_\eta(d, X_i))_{i=1}^n$.
\STATE (b)
Draw Bayesian bootstrap weights $M^s_{ni}=e^s_i/\sum_{j=1}^n e^s_j$ where $e_i^s \stackrel{iid}{\sim} \textup{Exp}(1)$, $i=1,\dots,n$.
\STATE (c) Calculate the corrected posterior draw for the ATE:
\begin{equation}\label{recentered_bay_est}
\check{\tau}_\eta^{s}=\tau_\eta^{s}-\widehat{b}^s_{\eta},
\end{equation}
\begin{equation}\label{debiased_bay_est}
\tau_\eta^{s}= \sum_{i=1}^n M^s_{ni}\big(m^s_\eta(1,X_i)-m^s_\eta(0,X_i)\big)\quad \text{and}\quad \widehat{b}^s_{\eta}=\frac{1}{n}\sum_{i=1}^n \boldsymbol{\tau}[m_\eta^s-\widehat m](Z_i),
\end{equation}
where $\boldsymbol{\tau}[m](z):=m(1,x)-m(0,x)+\widehat{\gamma}(d,x)(y-m(d,x))$.
\ENDFOR
\STATE \textbf{Output: $\{\check{\tau}_\eta^{s}:s=1,\ldots,S\}$}
\end{algorithmic}
\end{algorithm}
Given the draws from the corrected posterior calculated in Algorithm \ref{algorithm}, 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
\begin{equation*}
\mathcal{C}_n(\alpha)=\big\{\tau: q_n(\alpha/2)\leq \tau \leq q_n(1-\alpha/2)\big\},
\end{equation*}
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 \eqref{Riesz_est}, 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 \cite{ghosal2006binary}. Section \ref{sec:implement:GP} 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{appendix:simu}.
\begin{remark} [Bayesian bootstrap]\label{rem:BB} Under unconfoundedness and the reparametrization in \eqref{repar}, the ATE can be written as $\tau_{\eta}=\int [ \Psi\left(\eta^m(1,x)\right)- \Psi\left(\eta^m(0,x)\right)]\,\mathrm{d}F_\eta(x)$.
With independent priors on $\eta^m$ and $F_\eta$, their posteriors also become independent. It is thus sufficient to consider the posterior for $\eta^m$ and $F_\eta$ separately. We place a Dirichlet process prior for $F_\eta$ with the base measure to be zero. Consequently, the posterior law of $F_\eta$ coincides with the Bayesian bootstrap \citep{rubin1981bayesian}; also see \cite{chamberlain2003bayesian}. One key advantage of the Bayesian bootstrap is that it allows us to incorporate a broad class of data generating processes, whose posterior can be easily sampled. Replacing $F_\eta$ by the standard empirical cumulative distribution function does not provide sufficient randomization of $F_\eta$, as it yields an underestimation of the asymptotic variance;
see \cite[p. 3008]{ray2020causal}. In principle, one could consider other types of bootstrap weights; however, these generally do not correspond to the posterior of any given prior distribution.
\end{remark}
\section{Main Theoretical Results}\label{sec:asympt}
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:
\begin{equation}\label{def:est:chi}
\widehat{\tau}=\tau_0+\frac{1}{n}\sum_{i=1}^n \widetilde{\tau}_0(Z_i)+o_{P_0}(n^{-1/2}),
\end{equation}
where $\widetilde{\tau}_0=\widetilde{\tau}_{\eta_0}$ is the efficient influence function in accordance with \eqref{eif_ate}. 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 \cite{hahn1998role}.
\subsection{Least Favorable Direction}
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
\begin{eqnarray}\label{submodel}
\pi_t(\cdot ) = \Psi(\eta^{\pi}+t\mathfrak{p})(\cdot),~~
m_t(\cdot) = \Psi(\eta^m+t\mathfrak{m})(\cdot ),~~
f_t(\cdot )= \frac{f(\cdot)e^{t\mathfrak{f}(\cdot)}}{\int e^{t\mathfrak{f}(x)}f(x)\,\mathrm{d}x},
\end{eqnarray}
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 \emph{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:
\begin{equation}\label{score_ate}
B_{\eta}(\mathfrak{p},\mathfrak{m},\mathfrak{f})(z)= B_{\eta}^{\pi}\mathfrak{p}(z) + B_{\eta}^{m}\mathfrak{m}(z) + B_{\eta}^{f}\mathfrak{f}(z),
\end{equation}
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 \citet[p.370]{ghosal2017fundamentals}. We immediately obtain the following.
\begin{lemma}\label{lemma:lfd}
Consider the submodel \eqref{submodel}.
Let Assumption \ref{Ass:unconfounded} hold for $P_\eta$ with any $\eta$ under consideration, then
the least favorable direction for estimating the ATE parameter in \eqref{ate} is:
\begin{equation}\label{lfd}
\xi_{\eta}(d,x)= \left(0,\gamma_\eta(d,x), m_\eta(1,x)-m_\eta(0,x)-\tau_\eta\right),
\end{equation}
where the Riesz representer $\gamma_\eta$ is given in \eqref{riesz:def}.
\end{lemma}
Lemma \ref{lemma:lfd} motivates the adjustment of the prior distribution as considered in our Bayesian procedure in Section \ref{sec:method_outline}. 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{lemma:lfd} extends the result in Section 2.1 in \cite{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.
\subsection{Assumptions for Inference}
We now provide additional notations and assumptions. The posterior distribution plays an important role in the following analysis and is given by
\begin{equation*}
\Pi\left((\pi,m)\in A, F\in B|Z^{(n)}\right)=\int_{B}\frac{\int_{A}\prod_{i=1}^{n}p_{\pi,m}(Y_i, D_i |X_i) \,\mathrm{d}\Pi(\pi,m)}{\int \prod_{i=1}^{n}p_{\pi,m}(Y_i, D_i | X_i) \,\mathrm{d}\Pi(\pi,m)}\mathrm{d}\Pi(F|X^{(n)})
\end{equation*}
where $p_{\pi,m}$ denotes the conditional density of $(Y_i, D_i) $ given $X_i$, given by \eqref{x_den} 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{x_den}) 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$.
\begin{assumption}\label{Assump:Rate}[Rates of Convergence]
The estimators $\widehat \pi$ and $\widehat m$, which are based on an auxiliary sample independent of $Z^{(n)}$, satisfy $\Vert \widehat{\pi}-\pi_0\Vert_{2, F_0}=O_{P_0}(r_n) $ and for $d\in\{0,1\}$:
\begin{equation*}
\Vert \widehat{m}(d,\cdot)-m_0(d,\cdot)\Vert_{ 2,F_0}=O_{P_0}(\varepsilon_n)~\text{ and }~\sup_{\eta\in\mathcal{H}_n}\Vert m_\eta(d,\cdot)-m_0(d,\cdot)\Vert_{2,F_0}\lesssim \varepsilon_n,
\end{equation*}
where $\max\{\varepsilon_n, r_n\}\to 0$ and $\sqrt{n}\,\varepsilon_nr_n\to 0$. Further, $\Vert \widehat\gamma\Vert_{\infty}=O_{P_0}(1)$.
\end{assumption}
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)$.
\begin{assumption}\label{Assump:Donsker}[Complexity]
For $\mathcal{G}_n=\{\bar{m}_{\eta}(\cdot):\eta\in \mathcal{H}_n\}$ it holds
$\sup _{\bar{m}_{\eta} \in \mathcal{G}_n}\left|(\mathbb{P}_n-P_0) \bar{m}_{\eta}\right| =o_{P_0}(1)$
and
\begin{align}\label{NewSE}
\sup_{\eta\in\mathcal{H}_n}\left|\mathbb{G}_n\left[\left(\widehat\gamma-\gamma_0\right)(m_\eta-m_0)\right]\right|&=o_{P_0}(1).
\end{align}
\end{assumption}
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$.
\begin{assumption}\label{Assump:Prior}[Prior Stability]
For $d\in\{0,1\}$, $W^m(d,\cdot)$ is a continuous stochastic process independent of the normal random variable $\lambda\sim N(0,\sigma_n^2)$, where $\sigma_n\lesssim 1$, $n\sigma^2_{n}\to\infty$ and that satisfies: (i)
$\Pi\left(\lambda:|\lambda|\leq u_n\sigma_n^2\sqrt{n}\mid Z^{(n)}\right)\to_{P_0}1$,
for some deterministic sequence $u_n\to 0$ and (ii)
$\Pi\left((w,\lambda):w+(\lambda+tn^{-1/2})\widehat{\gamma}\in\mathcal{H}_n^m\mid Z^{(n)} \right)\to_{P_0}1$
for any $t\in\mathbb R$.
\end{assumption}
\textit{Discussion of Assumptions:}
Assumption~\ref{Assump:Rate} 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
\cite{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 \cite{ray2020causal}. The rate conditions in Assumption 2 also resemble conditions (i) and (ii) of Theorem~1 of \cite{farrell2015} in the context of frequentist
estimation. Remark~\ref{rem:DR} illustrates that under classical smoothness assumptions, this assumption is less restrictive than the method of \cite{ray2020causal} or
other approaches for semiparametric estimation of ATEs as found in
\cite{chen2008semiparametric} or \cite{farrell2021deep}. Assumption~\ref{Assump:Prior} incorporates Conditions (3.9) and (3.10) from Theorem~2 in \cite{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{sec:gauss:prior} for concreteness.
Assumption~\ref{Assump:Donsker} restricts the functional class
$\mathcal{G}_{n}$ to form a $P_{0}$-Glivenko--Cantelli class; see Section~2.4 of \cite{van1996empirical}.
This imposes a new stochastic equicontinuity condition, as (\ref{NewSE}) restricts a product structure involving $\widehat\gamma$ and $m_{\eta}$, which further relaxes the corresponding condition from \cite{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{rem:DR}). 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{Assump:Donsker} is also used in the frequentist
literature; see Section~2 of Benkeser, Carone, van der Laan, and
Gilbert (\citeyear{benkeser2017doubly}). Nonetheless, the
technical argument differs substantially from the frequentist's study,
because we mainly need the condition (\ref{NewSE}) 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.
\subsection{A Double Robust Bernstein-von Mises Theorem}
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
\begin{equation}
d_{BL}(P,Q)=\sup_{f\in BL(1)}\left| \int_{\mathcal{Z}}f(\mathrm{d}P-\mathrm{d}Q)\right|,
\end{equation}
where
\begin{equation*}
BL(1)=\left\{f:\mathcal{Z}\mapsto\mathbb{R}, \sup_{z\in\mathcal{Z}}|f(z)|+\sup_{z\neq z'}\frac{|f(z)-f(z')|}{\|z-z'\|_{\ell_2}}\leq 1 \right\}.
\end{equation*}
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{sec:method_outline}, to approximate this posterior distribution, as well as the resulting point estimator and credible set.
\begin{theorem}\label{thm:BvM}
Let Assumptions \ref{Ass:unconfounded}--\ref{Assump:Prior} hold. Then we have
\begin{equation*}
d_{BL}\left(\mathcal{L}_{\Pi}(\sqrt{n}(\tau_\eta-\widehat{\tau}-b_{0,\eta})|Z^{(n)}), N(0,\textsc v_0) \right)\to_{P_0} 0,
\end{equation*}
where $b_{0,\eta}:= \mathbb{P}_n[\gamma_0(m_0-m_{\eta})-(\bar{m}_0-\bar{m}_{\eta})] $.
\end{theorem}
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 \cite{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 \eqref{debiased_bay_est}, results in a feasible Bayesian procedure that satisfies the BvM theorem under double robustness, as demonstrated below.
\begin{theorem}\label{thm:Debias}
Let Assumptions \ref{Ass:unconfounded}--\ref{Assump:Prior} hold. Then we have
\begin{equation*}
d_{BL}\left(\mathcal{L}_{\Pi}(\sqrt{n}(\tau_\eta-\widehat{\tau}-\widehat{b}_\eta)|Z^{(n)}), N(0,\textsc v_0) \right)\to_{P_0} 0.
\end{equation*}
\end{theorem}
We now show how Theorem \ref{thm:Debias} 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$.
\begin{corollary}\label{cor:CI}
Let Assumptions \ref{Ass:unconfounded}--\ref{Assump:Prior} hold. Then under $P_0$, we have
\begin{equation}
\sqrt{n}\left(\overline{\tau}_{\eta}-\tau_0\right)\Rightarrow N(0,\textsc v_0).
\end{equation}
Also, for any $\alpha\in(0,1)$ we have
$P_0\big(\tau_0\in \mathcal{C}_n(\alpha)\big) \to 1-\alpha$.
\end{corollary}
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 \cite{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 \cite{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{Assump:Donsker}, 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.
\begin{remark}[Connection with frequentist robust estimation]
In our BvM theorem, we do not restrict the centering estimator $\widehat{\tau}$, as long as it admits the linear representation as in \eqref{def:est:chi}.
A popular frequentist estimator for the ATE that achieves double robustness is
\begin{equation}\label{freq_DR}
\widehat{\tau}= n^{-1}\sum_{i=1}^n \big(\widehat m(1,X_i)-\widehat m(0,X_i)\big) + n^{-1}\sum_{i=1}^n\widehat\gamma(D_i,X_i)\big(Y_i-\widehat m(D_i,X_i)\big)
\end{equation}
based on frequentist-type pilot estimators $\widehat{m}$ of the conditional mean function $m_0$ and $\widehat{\gamma}$ of the Riesz representer $\gamma_0$; see \cite{robins1995} and more recently \cite{ChernozhukovNeweySingh2020a,ChernozhukovNeweySingh2020b}.
The double robust or double machine learning estimator \eqref{freq_DR} recenters the plug-in type functional by an explicit correction factor that depends on the Riesz representer.\footnote{Another popular method in the statistics literature is the targeted learning approach \citep{vanderLaan2011tl,benkeser2017doubly}. }
Our main result establishes the asymptotic equivalence of our estimator and \eqref{freq_DR}. This not only offers frequentist validity to our Bayesian procedure but also provides a Bayesian interpretation for doubly robust frequentist methods.
\end{remark}
\begin{remark}[Parametric Bayesian Methods]
A couple of recent papers propose doubly robust Bayesian recipes for ATE inference, under parametric model restrictions. \cite{saarela2016double} considered a Bayesian procedure based on an analog of the double robust frequentist estimator given in Equation \eqref{freq_DR}, replacing the empirical measure with the Bayesian bootstrap measure. However, there was no formal BvM theorem presented therein.
Another recent paper by \cite{yiu2020unequal} explored Bayesian exponentially tilted empirical likelihood with a set of moment constraints that are of a double-robust type. They proved a BvM theorem for the posterior constructed from the resulting exponentially tilted empirical likelihood under parametric specifications.
\cite{luo2023semiparametric} provided Bayesian results for ATE estimation in a partial linear model, which implies homogeneous treatment effects. They also assign parametric priors to the propensity score. Their BvM Theorem allows for misspecification only in a parametric nonlinear component of the outcome equation. It is not clear how to extend their analysis to incorporate flexible nonparametric modeling strategies.
\end{remark}
\section{Illustration using Squared Exponential Process Priors}\label{sec:gauss:prior}
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 \citep{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.
\subsection{Asymptotic Results under Primitive Conditions}
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 \cite{rassmusen2006gaussian} and \cite{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 \cite{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
\begin{equation}\label{RescaleRate}
a_n\sim n^{1/(2s_m+p)}(\log n)^{-(1+p)/(2s_m+p)}.
\end{equation}
Under (\ref{RescaleRate}), 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 \cite{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.
\begin{proposition}\label{prop:exponential}
Let Assumption \ref{Ass:unconfounded} hold.
The estimator $\widehat\gamma$ satisfies $\|\widehat\gamma\|_\infty=O_{P_0}(1)$ and $\|\widehat{\gamma}-\gamma_0\Vert_\infty= O_{P_0}\big((n/\log n)^{-s_\pi/(2s_\pi+p)}\big)$ for some $s_\pi>0$. Suppose $m_0(d,\cdot)\in \mathcal{C}^{s_m}([0,1]^{p})$ for $d\in\{0,1\}$ and some $s_m>0$ with $\sqrt{s_\pi \, s_m}>p/2$. Also, $\|\widehat{m}(d,\cdot)-m_0(d,\cdot)\Vert_{2, F_0}= O_{P_0}\big((n/\log n)^{-s_m/(2s_m+p)}\big)$.
Consider the propensity score-dependent prior on $m$ given by $m(d,x) = \Psi\left(W^m(d,x) + \lambda\,\widehat \gamma(d,x)\right)$, where $W^m(d,\cdot)$ is the rescaled squared exponential process for $d\in\{0,1\}$, with its rescaling parameter $a_n$ of the order in \eqref{RescaleRate} and
$\left(n/\log n\right)^{-s_m/(2s_m+p)}\lesssim u_n\sigma_n$ for some deterministic sequence $u_n\to 0$, and $\sigma_n\lesssim 1$.
Then, the corrected posterior distribution for the ATE satisfies Theorem \ref{thm:BvM}.
\end{proposition}
\begin{remark}[Double Robust H\"older Smoothness]\label{rem:DR}
Proposition \ref{prop:exponential} requires $\sqrt{s_\pi \, s_m}>p/2$, which represents a trade-off between the smoothness requirement for $m_0$ and $\pi_0$. This encapsulate the \textit{double robustness}; i.e., a lack of smoothness of the conditional mean function $m_0$ can be mitigated by exploiting the regularity of the propensity score and vice versa. Referring to the H\"older class $\mathcal{C}^{s_m}([0,1]^{p})$, its complexity measured by the bracketing entropy of size $\varepsilon$ is of order $\varepsilon^{-2\upsilon}$ for $\upsilon=p/(2s_m)$. One can show that the key stochastic equicontinuity assumption in \cite{ray2020causal}, that is, their condition (3.5), is violated by exploring the Sudkov lower bound \citep{han2021set} when $\upsilon>1$ or equivalently when $s_m<p/2$. In contrast, our framework accommodates this non-Donsker regime as long as $\sqrt{s_\pi \, s_m}>p/2$, which enables us to exploit the product structure and a fast convergence rate for estimating the propensity score. Our methodology is not restricted to the case where propensity score belongs to a H\"older class per se. For instance, under a parametric restriction (such as in logistic regression) or an additive model with unknown link function, the possible range of the posterior contraction rate $\varepsilon_n$ for the conditional mean function can be substantially enlarged.
In the case $ s_m>p/2$, the bias term becomes asymptotically negligible, i.e., $b_{0,\eta}=o_{P_0}(n^{-1/2})$. This allows for smoothness robustness only with respect to the propensity score and is also known as single robustness. In this case, no posterior correction is required, see \cite{ray2020causal}.
\end{remark}
\subsection{Implementation Details}\label{sec:implement:GP}
We provide details on the Gaussian process prior placed on $\eta^m(d,x)$ and its posterior computation. Algorithm \ref{algorithm} 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 \citep[p.83]{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 \cite{ray2019debiased} and \cite{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{algorithm}, satisfies the conditions $\sigma_n\lesssim 1$ and $n\sigma^2_{n}\to\infty$ in Assumption \ref{Assump:Prior}, with probability approaching one. It is similar to the choice suggested by \citet[page 6]{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{appendix:simu} 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{appendix:LPA}.
For the implementation of the pilot estimator $\widehat\gamma$ given in \eqref{Riesz_est}, we recommend Logistic Lasso for the propensity score, with the penalty parameter chosen by cross-validation \citep{friedman2010regularization}.
As a pilot estimator $\widehat{m}$ in Algorithm \ref{algorithm} 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{algorithm}, 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{prop:exponential}, 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 \cite{ghosal2017fundamentals}.
\section{Numerical Results}
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 \cite{lalonde1986}, \cite{dehejia1999causal}, \cite{abadie2011bias}, and \cite{armstrong2021finite}, among others. We refer readers to \cite{lalonde1986}, and \cite{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}.}
\subsection{Simulations}\label{sec:simu}
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 \cite{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 \cite{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{algorithm}, using $S=5000$ posterior draws and the pilot estimator $\widehat \gamma$ and $\widehat m$, as detailed at the end of Section \ref{sec:implement:GP}.
We compare DR Bayes to two other Bayesian procedures:
First, we consider the prior adjusted Bayesian method (PA Bayes) proposed by \cite{ray2020causal}, which
constructs the point estimate and credible interval based on $\tau_{\eta}^s$ in (\ref{debiased_bay_est}). 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 \cite{abadie2011bias}, which adjusts for differences in covariate values through regression.
DR TMLE corresponds to the doubly robust targeted maximum likelihood estimator by \cite{benkeser2017doubly}.
DML refers to the double/debiased machine learning estimator from \cite{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 \cite{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{\cite{crump2009dealing} suggested a simple rule of thumb with a threshold of $t=0.10$, while \cite{athey2021using} used $t=0.05$. Applying the optimal trimming rule proposed by \cite{crump2009dealing} to our simulated samples yields an average optimal trimming threshold $0.073$.}
\begin{table}[H]
\centering
\caption{Simulation results using WGAN-generated data: trimming is based on the estimated propensity score within $[t,1-t]$, $\bar n =$ the average sample size after trimming, CP = coverage probability of $95\%$ credible/confidence interval, CIL = average length of the $95\%$ credible/confidence interval.
\qquad}\label{tab:simu_1}
\vskip.15cm
{\footnotesize
\begin{tabular}{lccccccccccccc}\toprule
\multicolumn{1}{c}{Methods}&\multicolumn{1}{c}{}& \multicolumn{1}{c}{Bias}&\multicolumn{1}{c}{ CP}&\multicolumn{1}{c}{CIL}&\multicolumn{1}{c}{}&\multicolumn{1}{c}{Bias}& \multicolumn{1}{c}{CP}&\multicolumn{1}{c}{CIL}&\multicolumn{1}{c}{}&\multicolumn{1}{c}{Bias}& \multicolumn{1}{c}{CP}&\multicolumn{1}{c}{CIL}\\
\cline{1-1}\cline{3-5}\cline{7-9}\cline{11-13}
\multicolumn{1}{c}{}&\multicolumn{1}{c}{}&\multicolumn{3}{c}{$t = 0.10 (\bar n = 240)$ }&\multicolumn{1}{c}{}& \multicolumn{3}{c}{$t =0.05 (\bar n =363)$}&\multicolumn{1}{c}{}&\multicolumn{3}{c}{$t = 0.01 (\bar n =664)$} \\
\cline{3-13}
Bayes&& -0.040 & 0.683 & 0.147&&-0.010 & 0.841 &0.149&& -0.006 & 0.911 & 0.120 \\
PA Bayes && -0.008 & 0.981 & 0.260 && 0.033 & 0.949 & 0.254 &&0.047 & 0.897 & 0.308 \\
DR Bayes && -0.024 & 0.983 & 0.223 && 0.014 & 0.970 & 0.221 &&0.023 & 0.952 & 0.258 \\
\cline{1-5}\cline{7-9}\cline{11-13}
Match && 0.027 & 0.933 & 0.334 && 0.048 &0.908 & 0.323&& 0.033 & 0.965 & 0.323 \\
Match BC && 0.040 &0.880 & 0.347 && 0.065 & 0.816 & 0.334&& 0.083 & 0.804 & 0.339 \\
DR TMLE&& 0.015 & 0.832 & 0.300 && 0.039 & 0.746 & 0.282&&0.039 & 0.668 & 0.242 \\
DML && 0.045 & 0.927 & 0.524&& 0.052 & 0.870 & 0.393&&0.054 & 0.918 & 0.522 \\
\bottomrule
\end{tabular}}
\end{table}
Table \ref{tab:simu_1} 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{appendix:simu} 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{tab:simu_1} 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 \cite{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.
\subsection{An Empirical Illustration}
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 \cite{crump2009dealing} yields an optimal threshold of $0.064$.} The ATE point estimates and confidence intervals are presented in Table \ref{tab:PSID_1}.
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]$.
\begin{table}[H]
\centering
\caption{Estimates of ATE for the Lalonde--Dehejia--Wahba data: trimming is based on the
estimated propensity score within $[t,1-t]$, $\bar n =$ sample size after trimming. ATE = point estimate, $95\%$ CI = $95\%$ credible/confidence interval, CIL = $95\%$ credible/confidence interval length.
\qquad}\label{tab:PSID_1}
\vskip.15cm
{\footnotesize
\addtolength{\tabcolsep}{-1.5pt}
\begin{tabular}{lccc|ccc|ccc}\toprule
\multicolumn{1}{c}{Methods }&\multicolumn{3}{c}{$t = 0.10 (\bar n=245) $}&\multicolumn{3}{c}{$t = 0.05 (\bar n=398) $}&\multicolumn{3}{c}{ $t =0.01 (\bar n=740)$}\\
\cline{2-10}
\multicolumn{1}{c}{}& \multicolumn{1}{c}{ATE}&\multicolumn{1}{c}{$95\%$ CI}&\multicolumn{1}{c}{CIL}&\multicolumn{1}{c}{ATE}& \multicolumn{1}{c}{$95\%$ CI}&\multicolumn{1}{c}{CIL}&\multicolumn{1}{c}{ATE}& \multicolumn{1}{c}{$95\%$ CI}&\multicolumn{1}{c}{CIL}\\
\cline{1-10}
Bayes& 0.213 & [0.120, 0.301] & 0.181 & 0.214 & [0.132, 0.292] & 0.161 &0.198 & [0.140, 0.251] & 0.112 \\
PA Bayes& 0.158 & [0.019, 0.288] & 0.270 & 0.170 & [0.045, 0.281] & 0.236 & 0.090 & [-0.078, 0.233] & 0.311\\
DR Bayes & 0.178 & [0.061, 0.293] & 0.231& 0.184 & [0.064, 0.294] & 0.230& 0.121 & [-0.031, 0.250] & 0.281\\
\cline{1-10}
Match & 0.188 & [0.022, 0.355] & 0.333 & 0.140 & [-0.029, 0.309] & 0.338 & 0.079 & [-0.111, 0.269] & 0.380\\
Match BC& 0.157 & [-0.006, 0.321] & 0.327 & 0.145 & [-0.021, 0.310] & 0.331 & 0.180 & [-0.004, 0.365] & 0.369\\
DR TMLE& -0.023 & [-0.171, 0.125] & 0.296 & 0.073 & [-0.074, 0.220] & 0.294 & 0.071 & [-0.146, 0.289] & 0.435\\
DML& 0.172 & [0.018, 0.327] & 0.308 &0.150 & [-0.010, 0.310] & 0.320 & 0.258 & [-0.183, 0.699] & 0.882\\
\bottomrule
\end{tabular}}
\end{table}
As we see from Table \ref{tab:PSID_1}, 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.
\section{Extensions}\label{sec:extension}
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{appendix:proof_exponent} and \ref{appendix:exponent}. Additionally, we outline extensions to other causal parameters of interest.
\subsection{A Single-parameter Exponential Family}\label{sec:single}
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
\begin{equation}\label{condpdf}
f_{Y|D,X}(y\mid d,x) = c(y)\exp\left[q(m(d,x))ay-A(m(d,x))\right],
\end{equation}
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{condpdf}) 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 \eqref{condpdf} 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
\begin{equation}\label{x_den_exp}
p_{\pi,m,f}(y,d,x)=\pi(x)^d (1-\pi(x))^{1-d}c(y)\exp\left[q(m(d,x))ay-A(m(d,x))\right]f(x).
\end{equation}
We consider the same reparametrization of $(\pi, m, f)$ as in \eqref{repar} 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.
\begin{lemma}\label{lemma:lfd_exp}
Let Assumption \ref{Ass:unconfounded} hold for $P_{\eta}$ for any $\eta$ under consideration. Then, for the joint distribution \eqref{x_den_exp} and the submodel $t\mapsto \eta_t$ defined by the path $m_t(d,x) = q^{-1}(\eta^m+t\mathfrak{m})(d,x)$ with $(\pi_t, f_t)$ as defined in \eqref{submodel},
the least favorable direction for estimating the ATE parameter in \eqref{ate} is:
\begin{equation}\label{lfd_exp}
\xi_{\eta}(d,x)= \left(0,\frac{1}{a}\gamma_\eta(d,x), m_\eta(1,x)-m_\eta(0,x)-\tau_\eta\right),
\end{equation}
where the Riesz representer $\gamma_\eta$ is given in \eqref{riesz:def}.
\end{lemma}
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{lemma:lfd}.
To implement the double robust Bayesian procedure for general outcomes, one can still follow Algorithm \ref{algorithm}, 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{algorithm} 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 \cite{ray2019debiased}.
Proposition \ref{prop:OnExponential} in the Supplemental Material provides primitive conditions for the BvM Theorem to hold under double robust smoothness conditions.
\subsection{Multinomial Outcomes}\label{sec:multi}
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:
\begin{align*}
\Psi_0\left(\eta^{m_1},\cdots,\eta^{m_J}\right)=\frac{1}{1+\sum_{l=1}^J \exp(\eta^{m_l})}
\quad \text{and}\quad
\Psi_j\left(\eta^{m_1},\cdots,\eta^{m_J}\right)=\frac{\exp(\eta^{m_j})}{1+\sum_{l=1}^J \exp(\eta^{m_l})},
\end{align*}
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.
\begin{lemma}\label{lemma:lfd:multinomial}
Consider the submodel $t\mapsto \eta_t$ defined by the path $ m_{t, j}(d,x) = \Psi(\eta^{m_j}+t\mathfrak{m}_j)(d,x)$, $1\leq j\leq J$, with $(\pi_t, f_t)$ as defined in \eqref{submodel}.
Let Assumption \ref{Ass:unconfounded} hold for $P_{\eta}$ for any $\eta$ under consideration, then the least favorable direction for estimating the ATE parameter is:
\begin{equation*}
\xi_{\eta}(d,x)= \left(0,\gamma_\eta(d,x), 2\gamma_\eta(d,x), \ldots, J\gamma_\eta(d,x), m_\eta(1,x)-m_\eta(0,x)-\tau_\eta\right),
\end{equation*}
where the Riesz representer $\gamma_\eta$ is given in \eqref{riesz:def}.
\end{lemma}
We emphasize that the least favorable direction calculation is not a trivial extension of \cite{hahn1998role} or \cite{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 \cite{hahn1998role}. In the proof of Lemma \ref{lemma:lfd:multinomial}, 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{lemma:lfd:multinomial} 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$:
\begin{align*}
m_{\eta, j}(d,x) = \Psi_j\left(\eta^{m_1},\cdots,\eta^{m_J}\right)(d,x)\qquad \text{and}\qquad \eta^{m_j}(d,x)=W^{m_j}(d, x) + \lambda\,j \widehat \gamma(d,x),
\end{align*}
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{sec:method_outline} using $m_\eta(d,x)=\sum_{j=0}^J j\, m_{\eta,j}(d,x)$.
\subsection{Other Causal Parameters}\label{sec:other_param}
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{appendix:lfd}.
\begin{example}[Average Policy Effects]
The policy effect from changing the distribution of $X$ is
$ \tau_\eta^{P}=\int m_{\eta}(x)\,\mathrm{d}(G_1(x)-G_0(x))$, where the known distribution functions $G_1$ and $G_0$ have their supports contained in the support of the marginal covariate distribution $F_{\eta}$.
Following the general setup, $ \psi(z,m_{\eta})=\psi(m_{\eta}):=\int m_{\eta}(x)\,\mathrm{d}(G_1(x)-G_0(x))$ with its Riesz representer
$\gamma_{\eta}^{P}(x)=(g_1(x)-g_0(x))/f_{\eta}(x)$, where $g_1$ and $g_0$ stand for the density function of $G_1$ and $G_0$, respectively. \end{example}
\begin{example}[Average Derivative]
For a continuous scalar (treatment) variable $D$, the average derivative is given by
$\tau_\eta^{AD}=\mathbb{E}_\eta\left[\partial_dm_\eta(D,X)\right]$, where $\partial_dm$ denotes the partial derivatives of $m$ with respect to the continuous treatment $D$. Thus, we have $\psi(Z,m_{\eta})=\partial_dm_{\eta}(D,X)$ with its Riesz representer given by $\gamma_\eta^{AD}(D,X)=\partial_d \pi_\eta(D,X)/\pi_\eta(D,X)$, where here $\pi_\eta$ denotes the conditional density function of $D$ given $X$.
\end{example}