Bayesian Semiparametric Causal Inference: Targeted Doubly Robust Estimation of Treatment Effects

Bayesian semiparametric causal inference: Targeted doubly robust estimation of treatment effects

\affil[1]{Department of Statistics, Texas A&M University} \makeatletter \makeatother

\footnotetext[1]{Corresponding author.} \footnotetext{{\it Email addresses:} \hyperlink{[email removed]}{[email removed]} (Gözde Sert), \hyperlink{[email removed]}{[email removed]} (Abhishek Chakrabortty), \hyperlink{[email removed]}{[email removed]} (Anirban Bhattacharya).}

{\bf Keywords:}Average treatment effect, Bayesian debiasing, hierarchical learning, high-dimensional nuisance, semiparametric Bayesian inference, summary statistics modeling.

Introduction

Inferring the causal effect of a treatment or exposure is central to many scientific disciplines. While randomized controlled trials are the gold standard for causal estimation, they are often infeasible due to ethical, logistical, or financial constraints. A common alternative is to use {\it observational} data, which is typically easier to obtain, but also requires careful methodology to handle potential confounding and high dimensionality issues, while ensuring robust (unbiased) estimation of causal estimands. Among these, the {\it average treatment effect} (ATE) is a key popular estimand, measuring the treatment's overall causal impact, and is widely adopted in various scientific disciplines.

Estimation of the ATE is naturally linked to semiparametric inference, as its identification involves infinite-dimensional nuisance parameters bang2005doubly. Most existing approaches are frequentist, such as propensity score adjustment or matching rosenbaum1983central, rosenbaum1984reducing, and doubly robust (DR) estimators robins1994estimation, robins1995semiparametric. Recently, {\it Bayesian} semiparametric methods for ATE estimation have gained attention ray2019debiased, ray2020semiparametric, hahn2020bayesian, antonelli2022causal, linero2022, luo2023semiparametric, breunig2025double. Traditional Bayesian methods marginalize out nuisance parameters to obtain a posterior of the target parameter and can achieve desirable contraction rates ghosal2017fundamentals. However, strong regularization often induces nuisance estimation bias bickel2012semiparametric, rivoirard2012bernstein, castillo2015bernstein that jeopardizes the validity of Bayesian inference for low-dimensional targets, such as the ATE.

Several strategies have been proposed to mitigate this bias. One line of research modifies or tailors priors to incorporate the propensity score and better align the prior with the semiparametric model structure ray2019debiased, ray2020semiparametric. Another applies posterior corrections or influence function hahn98 (IF)-based updates breunig2025double, yiu2025. A related method by antonelli2022causal constructs a posterior for the ATE by plugging nuisance posterior samples into the IF, followed by an additional variance correction for valid inference.

Building on recent advances, we propose the {\it doubly robust debiased Bayesian (DRDB) procedure}, which provides a principled and scalable solution to nuisance bias in high-dimensional or complex settings. DRDB departs from existing Bayesian methods in two key ways. First, it adopts a {\it targeted} modeling strategy that focuses on {\it summary statistics} informative about the ATE, rather than the full data distribution. Second, it introduces a {\it Bayesian debiasing} mechanism that {\it learns nuisance bias} directly from data, eliminating the need for prior modification or post hoc correction ray2020semiparametric, breunig2025double, yiu2025. By decoupling inference for the ATE from nuisance estimation, DRDB ensures robustness of the marginal posterior and serves as a {\it Bayesian analogue} of the frequentist double machine learning framework chernozhukov2018double, maintaining validity even under high-dimensional and/or misspecified models.

A prominent usage of summary statistics in Bayesian inference appears in the approximate Bayesian computation (ABC) literature to mitigate issues with a low acceptance rate drovandi2015. DRDB instead leverages them in a targeted manner to separate the ATE from nuisance bias. Its key component is a {\it retargeting step} that models the nuisance bias using weighted observables (an idea akin to importance sampling) which naturally incorporates the {\it propensity score} (PS) into the Bayesian framework. Although the PS plays a central role in the frequentist literature on DR estimation robins1994estimation, robins1995semiparametric, bang2005doubly, it has lacked a principled Bayesian counterpart li2023bayesian. DRDB fills this gap by {\it integrating the PS seamlessly} through its debiasing mechanism. Related work by sert2025 develops Bayesian inference via summary statistics for semi-supervised learning, which motivates the construction of DRDB. However, DRDB differs in two key respects: (i) DRDB uses a hierarchical model to learn the ATE directly, without relying on independence between data subsets, and (ii) it introduces a retargeting mechanism to identify and estimate nuisance bias appropriately.

Building on this bias estimation, DRDB integrates the bias into a {\it hierarchical} Bayesian framework: The posterior for the bias informs a conditional likelihood for the ATE, whose integration yields a valid marginal posterior for the ATE (see Equation (ref)). Another salient feature of DRDB is its use of {\it sample-splitting} and cross-fitting (CF) chernozhukov2018double. {\it Beyond} their traditional role in technical aspects, DRDB uses them as critical {\it methodological} tools to validate the debiasing step and {\it decouple} nuisance estimation from target inference. DRDB employs randomized splitting to obtain multiple subposteriors and aggregates them using a consensus Monte Carlo–type scheme scott2022bayes, producing a posterior that efficiently utilizes the {\it entire} data (see Section (ref)).

DRDB establishes a semiparametric {\it Bernstein-von Mises (BvM) result} for the marginal posterior of the ATE (Theorems (ref) and (ref)): When both nuisance models are well-specified and their posteriors contract at rates whose {\it product} is $o(n^{-1/2})$, the posterior concentrates around the true ATE at the parametric rate and is asymptotically Gaussian. In this case, the posterior mean is an asymptotically {\it efficient} estimator of the true ATE, converging at a $\sqrt{n}$-rate with asymptotic variance that achieves the semiparametric efficiency bound hahn98. Notably, the DRDB posterior depends on the nuisance posteriors only through their asymptotic limits, underscoring its robustness to nuisance modeling. Moreover, DRDB satisfies {\it Bayesian double robustness}: when only one nuisance model is well-specified (consistently estimated, while the other may be misspecified or slowly estimated), the {\it posterior remains consistent} for the ATE, contracting at the rate of the well-specified nuisance, extending the frequentist DR principle bang2005doubly to Bayesian inference (posteriors).

Finally, the key principles of DRDB (its debiasing mechanism, targeted modeling, and hierarchical learning strategy) extend naturally beyond the ATE, providing valid Bayesian inference for a broad class of causal estimands, including the average treatment effect on the treated (ATT), that on the control (ATC), and subgroup-specific effects. For clarity and brevity, the detailed extension of DRDB to these general estimands is presented in Section (ref) of the \hyperref[sec_supplementary]{Supplementary Material}.

The rest of this paper is organized as follows. Section (ref) introduces the basic setup and preliminaries. Section (ref) develops our proposed DRDB methodology, first for one counterfactual mean (Section (ref)), then for the ATE (Section (ref)). Section (ref) presents the technical details of the DRDB posterior, and the main theoretical results, including BvM results and Bayesian double robustness. Section (ref) reports finite-sample performance via simulations and data analysis. Section (ref) provides a concluding discussion. Extensions of our methodology, additional simulation results, and proofs and technical details are deferred to the \hyperref[sec_supplementary]{Supplementary Material} (Sections (ref)-(ref)).

Setup and preliminary causal assumptions

Data and notation

Let $T \in \{0, 1\}$ denote a binary treatment indicator; $Y \in \mathbb{R}$ denote the {\it observed} outcome, defined as: $Y = TY(1) + (1-T)Y(0)$, where $\{ Y(1), Y(0)\}$ are the {\it potential outcomes} rubin1974estimating, imbens2015causal under treatment ($T = 1$) and control ($T = 0$), respectively (i.e., $Y(t)$ is the outcome that would have been observed if $T = t$, possibly contrary to fact); and $\mathbf{X} \in \mathbb{R}^p$ denote the vector of covariates (or potential {\it confounders}). The {\it observed data} $\mathcal{D}$ consists of independent and identically distributed (i.i.d) observations $\mathbf{Z}_1, \dots, \mathbf{Z}_n$ of the random variable $\mathbf{Z} := (Y,\mathbf{X}, T)$ with support $\mathcal{Y} \times \mathcal{X} \times \{0, 1\}$ and underlying joint probability distribution (p.d.) $\mathbb{P}_{\mathbf{Z}}$. Also, the setting is throughout allowed to be (possibly) {\it high dimensional} (i.e., $p$ is allowed to grow with $n$).

Let $U$ be a random object with an underlying p.d. $\mathbb{P}_{U}$, and $f$ be a measurable $\mathbb{R}$-valued function of $U$. The expectation of $f(U)$ is defined as $\mathbb{E}_{U}\{f(U)\}:= \int f(u) d\mathbb{P}_{U}(u)$, whenever it exists. For any $d\geq 1$, $L_d(\mathbb{P}_{U})$ denotes the space of all $\mathbb{R}$-valued measurable functions of $U$ equipped with the norm $\|f \|_{L_d(\mathbb{P}_{U})}:= [\mathbb{E}_{U}\{f(U)^d\}]^{1/d}$. {\it We adopt the following Bayesian notation throughout: for a generic random object $\theta$, $\Pi_{\theta}$ denotes its posterior, $\underline{\theta}$ a posterior sample, and $\theta^\dagger$ its true value}.

Identification

The parameter of interest is the {\it average treatment effect (ATE), defined as:} $\Delta^\dagger:= \mu^\dagger(1) - \mu^\dagger(0)$, with $\mu^\dagger(t):= \mathbb{E}_{\mathbb{Z}}\{Y(t)\}$ for $t \in \{0, 1\}$, where the expectation is taken under the true p.d. $\mathbb{P}_{\mathbb{Z}}$ of $\mathbb{Z} := \{Y(1), Y(0), \mathbf{X}, T\}$. Since $\{Y(1), Y(0)\}$ cannot be jointly observed in the data, we impose standard causal assumptions to identify $\Delta^\dagger$ from the available data $\mathcal{D}$ rosenbaum1984reducing.

Assumption (ref)(a), known as {\it no unmeasured confounding (NUC)}, posits that the set of observed covariates captures all confounding factors affecting both treatment assignment and the potential outcomes. Assumption (ref)(b) imposes an {\it overlap} condition, ensuring that $\mathbf{X}$ in the treatment groups (i.e., $\mathbf{X} \,| \, T = t$) share sufficient common support for valid comparisons imbens2015causal.

\paragraph{Regression-based identification.} Define $m^\dagger_t(\mathbf{X}) \equiv m^\dagger(\mathbf{X}, t) := \mathbb{E}(Y(t) \mid \mathbf{X})$ as the {\it regression function} for treatment $t \in \{0,1\}$. Under Assumption (ref) (ignorability), it can be equivalently written as: $m_t^\dagger(\mathbf{X}) = \mathbb{E}(Y \mid \mathbf{X}, T = t)$. The ATE is then {\it identified} using the law of iterated expectations:

where we note that each $m_t^\dagger(\cdot)$ is {\it estimable} via a regression in the {\it observable} data on: $(Y, \mathbf{X}) \hspace{0.01in}| T = t$. Hence, the ATE $\Delta^\dagger$ is a {\it functional} of both $\mathbb{P}_{\mathbf{X}}$ and the {\it nuisance functions} $\overrightarrow{m}^\dagger \equiv \overrightarrow{m}^\dagger(\cdot)$, with $\overrightarrow{m}^\dagger:= (m_1^\dagger, m_0^\dagger)$, and this identification serves as a {\it foundation for our approach} to estimating $\Delta^\dagger$.

Methodology

Motivated by (ref), we propose a doubly robust debiased Bayesian (DRDB) procedure for estimating the ATE. To clarify the main steps of DRDB, we first present the methodology for a {\it single-arm:} $\mu_1^\dagger\equiv\mu^\dagger(1) = \mathbb{E}[Y(1)]$, and thereafter, extend it to the ATE in Section (ref). Importantly, estimating $\mu_1^\dagger$ is an interesting and non-trivial problem in its own right, as it corresponds to the mean of an outcome that is missing at random (MAR) within the {\it missing data} framework tsiatis2007semiparametric.

Let $\mathbb{K} \ge 2$ be a {\it fixed} integer. We randomly {\it split} $\mathcal{D}$ into $\mathbb{K}$ disjoint subsets $\{\mathcal{D}_k\}_{k=1}^\mathbb{K}$, each of equal size $n_{\mathbb{K}} := n/\mathbb{K}$, assuming without loss of generality that $n$ is divisible by $\mathbb{K}$. The corresponding index sets are denoted by $\{\mathcal{I}_k\}_{k=1}^\mathbb{K}$. For each $k \in \{1, \dots, \mathbb{K}\}$, define $\mathcal{D}_k^\- := \mathcal{D} \setminus \mathcal{D}_k$, which has size $n_{\mathbb{K}}^\- := n - n_{\mathbb{K}}$ and index set $\mathcal{I}_k^\-$. Let $(S, S^\-) := (\mathcal{D}_k, \mathcal{D}_k^\-)$ denote a generic pair of test and training datasets with corresponding index sets $(\mathcal{I}, \mathcal{I}^\-)$ for some $k \in \{1, \dots, \mathbb{K}\}$. For $t \in {0,1}$, let $(S_t, S_t^\-)$ denote the {\it subgroups} of $(S, S^\-)$ corresponding to treatment level $T = t$, so $S_1$ and $S_1^\-$ represent the respective {\it treated} subgroups, and $S_0$ and $S_0^\-$ represent the {\it control} subgroups. By construction, $S$ and $S^\-$ are {\it independent} ($S \!\perp \!\!\! \perp\! S^\-$), which is both {\it crucial} and {\it necessary} for the DRDB approach.

\paragraph{Motivating the DRDB procedure.} An intuitive approach to estimating $\mu_1^\dagger$, motivated by (ref), is to use a regression-based Bayesian ({\tt BREG}) procedure: Suppose the unknown nuisance function $m_1^\dagger(\cdot)$ is {\it learned from $S^\-$} via {\it any} suitable Bayesian regression method--parametric (like Bayesian ridge regression via Gaussian priors, or high dimensional sparse Bayesian linear regression using spike-and-slab type priors johnson2012bayesian) or nonparametric (such as Gaussian process regression williams1998prediction or Bayesian additive regression trees (BART) bart2010)--yielding a {\it posterior $\Pi_{m_1}$ for $m_1$}. For a sample $\underbar{m}_1 \sim \Pi_{m_1}$, one can treat $\{\underbar{m}_1(\mathbf{X}_i)\}_{i \in \mathcal{I}}$ as {\it derived i.i.d. samples in $S$}, targeting $\mu_1^\dagger$ through their mean. A standard Bayesian analysis, specifying a likelihood for this data and a prior on model parameters, then yields a posterior $\Pi_{\mathrm{reg}}$ for $\mu_1$.

Despite its intuitive appeal, BREG is highly {\it sensitive} to the quality of nuisance estimation: A misspecified nuisance model leads to an inconsistent posterior $\Pi_{\text{reg}}$ for $\mu_1$. Even with a correctly specified nuisance model, the posterior's {\it first-order} properties, such as its rate and shape, are {\it strongly} determined by the nuisance estimation bias: $\mathbb{E}_\mathbf{X}\{\underbar{m}_1(\mathbf{X}) - m_1^\dagger(\mathbf{X}) | \underbar{m}_1\}$. This makes the posterior overly dependent on the behavior of $\Pi_{m_1}$ and the choice of regression method, which in turn requires restrictive conditions to control the bias (e.g., in high dimensions) for achieving BvM-type results. These limitations motivate the key principles of our DRDB approach, which systematically eliminates this nuisance estimation bias within a Bayesian likelihood framework.

Doubly robust debiased Bayesian (DRDB) procedure for \texorpdfstring{$\mu^\dagger(1)$} {mu(1)}

DRDB is fundamentally a two-step approach. First, a Bayesian debiasing step learns and corrects for nuisance estimation bias within a Bayesian framework via a retargeting method. Second, a hierarchical learning framework learns the parameter of interest, $\mu_1^\dagger \equiv \mu^\dagger(1)$, after this bias has been addressed. We detail the steps in subsequent sections.

Adopting the notation from Section (ref), let $(S, S^\-)$ be a pair of test and training datasets. Assume the nuisance posterior $\Pi_{m_1} \equiv \Pi_{m_1}(\cdot; S^\-)$ for $m_1$ is obtained from $S^\-$ as before.

\paragraph{Debiasing step.} Let $\underbar{m}_1 \sim \Pi_{m_1} \equiv \Pi_{m_1}(\cdot; S^\-)$ be {\it one} sample independent (by design) of $S$. Using the regression-based representation of $\mu_1^\dagger$ in given (ref), we obtain the {\it debiased identification} of $\mu_1^\dagger$:

The term $b^\dagger(\underbar{m}_1)$ captures the nuisance estimation {\it bias} from replacing the true $m_1^\dagger$ with a random sample $\underbar{m}_1$. This bias is the primary source of the limitations of {\tt BREG} and serves as the central target of our {\it Bayesian debiasing} strategy. Its analysis and the validity of the debiased decomposition in (ref) crucially rely on the {\it independence condition} that ensures the distribution of $\mathbf{X} \in S$ in $\underbar{m}_1(\mathbf{X})$ is {\it unaffected} by that of $\underbar{m}_1 \sim {\Pi_{m_1}(\cdot;S^\-)}$ since $S^\- \!\perp \!\!\! \perp\! S$. To further analyze $b^\dagger(\underbar{m}_1)$, we write it as:

This formulation implies that if $Y(1)$ and $\mathbf{X}$ were observed for all units in $S$, one could directly estimate $b^\dagger(\underbar{m}_1)$ from $S$. However, both $\{Y(1), \mathbf{X}\}$ are {\it only} observed in the {\it treated subgroup:} $S_1$. Moreover, given $\underbar{m}_1 \sim \Pi_{m_1}$, the observables $\{Y - \underbar{m}_1(\mathbf{X})\} \in S_1$ target $\mathbb{E}_{(Y, \mathbf{X}) | T = 1}\{Y - \underbar{m}_1(\mathbf{X})\}$, rather than the desired bias $b^\dagger(\underbar{m}_1)$. To correct this discrepancy, a {\it retargeting} step is required in which the (derived) observations $\{Y - \underbar{m}_1(\mathbf{X})\} \in S_1$ are {\it reweighted} using a {\it density ratio} function. This adjustment ensures that the distribution of the weighted observations {\it aligns} with that of $(Y, \mathbf{X})$ in the whole population, rather than the conditional distribution given $T = 1$.

\paragraph{Retargeting bias via weighting.} Let $r_1^\dagger(\mathbf{X}) := f(\mathbf{X})/f_1(\mathbf{X} | T = 1)$ be the {\it density ratio} function, where $f(\cdot)$ is the density function (pdf) of $\mathbf{X}$ and $f_1(\cdot)$ is the conditional pdf given $T = 1$. Given $\underbar{m}_1 \sim \Pi_{m_1}$, we define the {\it weighted} observations $r_1^\dagger(\mathbf{X})\{Y - \underbar{m}_1(\mathbf{X})\}$ in $S_1$ and observe that:

This derivation shows that {\it unbiasedly} estimating the bias requires using the {\it weighted} observables $r_1^\dagger(\mathbf{X})\{Y - \underbar{m}_1(\mathbf{X})\}$ in $S_1$. It also clarifies that the bias $b^\dagger(\underbar{m}_1) \equiv b^\dagger(\underbar{m}_1, r_1^\dagger)$ should be viewed as a {\it functional} of {\it two} nuisances: $\underbar{m}_1$ and $r_1^\dagger \equiv r_1^\dagger(\cdot)$. To learn $b^\dagger(\underbar{m}_1, r_1^\dagger)$ from $S_1$, we must first estimate $r_1^\dagger$ by deriving a posterior from $S^\-$. This leads to the final analysis of the bias $b^\dagger(\underbar{m}_1, r_1^\dagger)$.

Let $\underbar{r}_1$ be {\it one} sample from the posterior $\Pi_{r_1} \equiv \Pi_{r_1}(\cdot; S^\-)$ of $r_1 \equiv r_1(\cdot)$, which we derive from a Bayesian {\it binary regression} method (e.g., Bayesian logistic regression or BART bart2010), as detailed in Remark (ref). Since $\underbar{r}_1$ is independent of $S$, substituting it into (ref) yields:

(ref) provides a full characterization of the bias. The second term in (ref): $\Gamma^\dagger(\underbar{m}_1, \underbar{r}_1) := \mathbb{E}[\{r_1^\dagger(\mathbf{X}) - \underbar{r}_1(\mathbf{X})\} \{m_1^\dagger(\mathbf{X}) - \underbar{m}_1(\mathbf{X})\} |\underbar{m}_1, \underbar{r}_1]$ is a {\it second-order} `bias of bias' (or `drift') term, arising from the {\it product} of the estimation errors for $r_1^\dagger$ and $m_1^\dagger$. We subsequently focus on {\it modeling and correcting the more tractable, first-order, bias} $b^\dagger(\underbar{m}_1, \underbar{r}_1)$. The second-order term $\Gamma^\dagger(\underbar{m}_1, \underbar{r}_1)$, while accounted for in the theoretical analysis of our eventual posterior, is not the primary debiasing target.

\paragraph{Targeted modeling strategy for bias.} Given $\underbar{m}_1 \sim \Pi_{m_1}$ and $\underbar{r}_1 \sim \Pi_{r_1}$, the bias $b^\dagger(\underbar{m}_1, \underbar{r}_1) = \mathbb{E}_{S_1}[\underbar{r}_1(\mathbf{X})\{Y - \underbar{m}_1(\mathbf{X})\} | \underbar{m}_1, \underbar{r}_1]$ can be viewed as a {\it functional} of the underlying distribution of $S_1$, specifically, relying on the {\it summary statistic} of the weighted observables $\underbar{r}_1(\mathbf{X})\{Y - \underbar{m}_1(\mathbf{X})\} \in S_1$. We can then construct a {\it working} likelihood based on these i.i.d. observables in $S_1$ and place a prior on the model parameters, yielding {\it a posterior $\Pi_{b_1}$ for $b_1 \equiv b(\underbar{m}_1, \underbar{r}_1)$}, as detailed in Proposition (ref).

A defining feature of the DRDB procedure, beyond its debiasing mechanism, is the {\it targeted} use of data. While traditional methods model the entire data ray2019debiased, ray2020semiparametric, breunig2025double, DRDB exclusively targets the parameters directly informative for $\mu_1^\dagger$. This targeted modeling strategy, combined with debiasing, forms the core of our DRDB approach, distinguishing it from existing Bayesian methodologies. Building on these, we next introduce the {\it hierarchical learning framework}, which facilitates a {\it construction} of a valid marginal posterior for $\mu_1^\dagger$, by leveraging in a novel way the conventional integral representation of the marginal posterior.

\paragraph{Hierarchical learning strategy.} If one had access to a joint posterior for $\{\mu_1, b(\underbar{m}_1)\}$, which would be the case in a traditional Bayesian framework, the marginal posterior for $\mu_1$ would be obtained by integrating out $b(\underbar{m}_1)$: $[\mu_1| S] = \int [\mu_1 | b(\underbar{m}_1), S] \hspace{0.05cm}[b(\underbar{m}_1) | S_1] \hspace{0.05cm}\mathrm{d} b(\underbar{m}_1)$. We instead take a fundamentally different perspective, using the right-hand side of this representation itself as the {\it defining principle} for constructing a posterior for $\mu_1$. In this formulation, $[\mu_1 | b(\underbar{m}_1), S]$ and $[b(\underbar{m}_1)| S]$ are individual building blocks which we aim to learn separately under our targeted modeling strategy, and mixing (integrating) over the uncertainty in $b(\underbar{m}_1)$ defines an idealized marginal posterior for $\mu_1$. Since the actual bias $b(\underbar{m}_1)$ involves the intractable second-order term in (ref), we achieve a key simplification by replacing the mixing variable $b(\underbar{m}_1)$ with its first-order {\it proxy} $b_1 \equiv b(\underbar{m}_1, \underbar{r}_1)$. Thus, we replace $[b(\underbar{m}_1)| S]$ with the posterior $\Pi_{b_1}$ for $b_1$ alluded to earlier, and construct a conditional posterior $\Pi_{\mu_1|b_1}$ for $\mu_1|b_1$; see discussion after Remark (ref). These switches allow us to mix over the proxy variable $b_1$, and define a valid probability measure

which serves as our {\it constructive definition} of a posterior for $\mu_1$. Here, $\Pi_{b_1}(\cdot; S)$ coincides with $\Pi_{b_1}(\cdot; S_1)$, as it is derived solely from $S_1$ under the targeted modeling strategy (see Proposition (ref)). This formulation is not only a practical construction of a marginal posterior but also represents a conceptually distinct perspective, motivated directly by the structure of the Bayesian integral, offering a new way to define posteriors when a full joint model is unavailable.

To implement this hierarchical construction, it {\it suffices to specify $\Pi_{\mu_1|b_1}$}. We formulate a {\it conditional likelihood} on $S$, which allows us to exploit the debiased representation in (ref) and adhere to our target-specific strategy. Specifically, we take a sample $\underbar{\it b}_1 \sim \Pi_{b_1}$ and model $\mu_1^\dagger - \underbar{b}_1$ conditional on $\underbar{b}_1$, using $S$. Given $\underbar{b}_1 \sim \Pi_{b_1}$, we have i.i.d. observables $\{\underbar{m}_1(\mathbf{X}_i)\} \in S$, which {\it target} $\mu_1^\dagger - \underbar{b}_1$ via their mean. Using these observables, we construct a working conditional likelihood with a corresponding prior, yielding a conditional posterior $\Pi_{\mu_1 |b_1}$ for $\mu_1|b_1$ (see Proposition (ref)). Finally, combining $\Pi_{\mu_1| b_1}$ and $\Pi_{b_1}$ via the integral (ref), we get a {\it marginal posterior $\Pi_{\mu_1} \equiv \Pi_{\mu_1}(\cdot; S)$ for $\mu_1$}.

The initial DRDB formulation relies on a single data split, $(S^\-, S)$, to ensure the independence required for our debiasing and targeted modeling strategy. The drawback, however, is a significant loss of efficiency from using only a fraction of the data for the final inference. To address this, we now detail a strategy to construct a final posterior for $\mu_1$ based on usage of the {\it full data} $\mathcal{D}$.

\paragraph{The final DRDB posterior with cross-fitting.} Leveraging the randomized sample-splitting presented in Section (ref), we can apply the DRDB procedure, as detailed in Steps (ref)--(ref), to {\it each} of the $\mathbb{K}$ test and training folds $\{(\mathcal{D}_k, \mathcal{D}_k^\-)\}_{k=1}^\mathbb{K}$. This yields {\it corresponding posteriors:} $\Pi_{\mu_1}^{(1)}, \ldots, \Pi_{\mu_1}^{(\mathbb{K})}$.

To efficiently use all available data, we {\it aggregate} these fold-specific posteriors into a {\it final posterior} for $\mu_1$ that incorporates information from all splits. Following the {\it consensus Monte Carlo (CMC)-type aggregation} strategy employed in sert2025, we define a new random variable, $\mu_1^{\mathrm{CF}}$, as the average of independent samples $\{\mu_1^{(k)}\}_{k = 1}^\mathbb{K}$ drawn from the respective posteriors $\{\Pi_{\mu_1}^{(k)}\}_{k=1}^\mathbb{K}$:

The resulting distribution, $\Pi_{\mu_1}^\mathrm{CF}$, serves as the {\it final DRDB posterior for $\mu_1$} and is a scaled convolution of the fold-specific posteriors $\{\Pi_{\mu_1}^{(k)}\}_{k=1}^\mathbb{K}$. This construction provides a principled and computationally efficient way to unify inference for $\mu_1^\dagger$ across all splits.

Although the combination step draws inspiration from CMC scott2022bayes, its goal here is quite different. We use sample splitting not for computational efficiency, but as a methodological necessity to create independent training and test sets that validate the debiased representation and enable targeted modeling. The aggregation strategy then provides a {\it Bayesian analogue of cross-fitting} (CF) chernozhukov2018double of {\it posteriors}, for semiparametric inference problems.

Generalized DRDB procedure for the ATE

Building on the DRDB procedure for $\mu_1^\dagger$, we now extend the method to our primary target, the ATE, $\Delta^\dagger$. Unlike frequentist ATE {\it point estimators}, which are simply the difference between the mean estimates for the two arms, Bayesian inference has {\it no} such direct analogue of `subtracting' {\it posterior distributions}. Hence, constructing a {\it valid} posterior for the ATE requires a more careful `first-principles' approach, accounting for the {\it full} posterior structure of the underlying components.

For clarity, we first detail the generalized DRDB procedure for one data split $(S, S^\-) = (\mathcal{D}_k, \mathcal{D}_k^\-)$. Second, we use the combination step from Section (ref) to aggregate the posteriors from all $\mathbb{K}$ folds.

\paragraph{Debiasing step for the ATE.} The regression-based identification in (ref) shows the ATE $\Delta^\dagger$ can be expressed as a {\it functional} of $\overrightarrow{m}^\dagger$ and $\mathbb{P}_\mathbf{X}$, denoted $\Delta^\dagger = \Delta^\dagger(\overrightarrow{m}^\dagger, \mathbb{P}_{\mathbf{X}})$. The nuisance function $\overrightarrow{m}^\dagger = (m_1^\dagger, m_0^\dagger)$ includes the unknown regression functions $m_1$ and $m_0$, both requiring estimation.

Define $m^\dagger(\cdot) := m_1^\dagger(\cdot) - m_0^\dagger(\cdot)$. Let $\underbar{m} \equiv \underbar{m}(\cdot)$ be {\it one} draw from the posterior $\Pi_{m} \equiv \Pi_{m}(\cdot; S^\-)$ obtained from $S^\-$, as in Remark (ref). Using the debiasing framework in Section (ref), we express $\Delta^\dagger$ in its {\it debiased} form as: $\Delta^\dagger = b^\dagger(\underbar{m}) + \mathbb{E}_{\mathbf{X} \in S}\{\underbar{m}(\mathbf{X}) | \underbar{m}\}$, where $b^\dagger(\underbar{m}):= \mathbb{E}_{\mathbf{X} \in S}\{m(\mathbf{X}) -\underbar{m}(\mathbf{X}) | \underbar{m}\}$, with the equality being valid due to the independence condition ($\underbar{m} \perp \!\!\! \perp S$). The term $b^\dagger(\underbar{m})$ is the expected bias from learning $m^\dagger(\cdot)$ and can be viewed as a function of $(b_1^\dagger, b_0^\dagger) \equiv (b^\dagger(\underbar{m}_1), b^\dagger(\underbar{m}_0))$, where $b^\dagger(\underbar{\it m}_t)$ is the bias for each arm $\mu^\dagger(t)$ for $t = 0,1$. To accurately model $b^\dagger(\underbar{m})$ and construct its posterior $\Pi_b \equiv \Pi_b(\cdot; S)$, we require a {\it joint learning} strategy: both biases, $b_1^\dagger$ and $b_0^\dagger$, must be learned together to produce a joint posterior for $(b_1, b_0)$, which defines a valid posterior for $b(\underbar{m})$.

\paragraph{Bias modeling for the ATE.} Following the bias analysis in Section (ref) and adopting the notational conventions introduced there for the {\it first-order bias} and the {\it drift term} (see Equations (ref)--(ref)), the bias $b^\dagger(\underbar{m}) \equiv b^\dagger(\underbar{m}, r^\dagger)$ can be decomposed into a first-order bias $b^\dagger(\underbar{m}, \underbar{r})$ and a drift term $\Gamma^\dagger(\underbar{m}, \underbar{r})$:

where $r^\dagger \equiv r^\dagger(\cdot) := (r_1^\dagger(\cdot), r_0^\dagger(\cdot))$, with $r_1^\dagger(\mathbf{X}) = p_1/e^\dagger(\mathbf{X})$ and $r_0^\dagger(\mathbf{X}) := (1-p_1)/\{1-e^\dagger(\mathbf{X})\}$, and with corresponding posterior draws $\underbar{r}_1:= \widehat{p}_1/\underbar{\it e} \sim \Pi_{r_1}$ and $\underbar{r}_0:= (1-\widehat{p}_1)/(1-\underbar{\it e}) \sim \Pi_{r_0}$ (see Remark (ref) for details), yielding a posterior sample $\underbar{r}:= (\underbar{r}_0, \underbar{r}_1)$ from $(\Pi_{r_0}, \Pi_{r_1})$. (ref) emphasizes that retargeting with the density ratio is {\it crucial} for accurate bias estimation. Similar to the one-arm case (see the discussions around (ref) and (ref)), we {\it focus} on modeling the {\it first-order} bias $b^\dagger(\underbar{m}, \underbar{r})$ in (ref) above, consistent with our main goal of debiasing. The second-order term $\Gamma^\dagger(\underbar{m}, \underbar{r})$, though not our primary debiasing target, is included in the theoretical analysis of our eventual posterior of $\Delta^\dagger$.

\paragraph{Posterior calculation for bias.} Recall that $b^\dagger(\underbar{m}, \underbar{r})$ can be expressed as function of $(b_1^\dagger, b_0^\dagger) \equiv (b^\dagger(\underbar{m}_1, \underbar{r}_1), b^\dagger(\underbar{m}_0, \underbar{r}_0))$, where $b^\dagger(\underbar{\it m}_t, \underbar{\it r}_t)$ denotes the {\it first-order} bias for each arm $t = 0, 1$, as defined in (ref). Thus, calculating the posterior $\Pi_b$ for $b$ reduces to obtaining the {\it joint} posterior $\Pi_{(b_1, b_0)}$ for $(b_1, b_0)$ from $S$. Notably, since the treated and control subsets, $S_1$ and $S_0$, are {\it independent} (by design), the {\it joint} posterior $\Pi_{(b_1, b_0)}$ can be factorized as the product of the {\it marginal} posteriors:

where $\Pi_{b_t}$ denotes the posterior of $b_t$ based on $S_t$ for $t = 0,1$. For explicit derivations, see Proposition (ref) in Section (ref). This factorization not only simplifies the analysis of $\Pi_b$, but also provides a straightforward sampling procedure: first, draw a sample $\underbar{b}_1$ from $\Pi_{b_1}$ and $\underbar{b}_0$ from $\Pi_{b_0}$, then define $\underbar{b} := \underbar{b}_1 - \underbar{b}_0$, yielding a {\it posterior sample from $\Pi_b$} for constructing the ATE posterior.

\paragraph{Hierarchical learning and posterior construction for the ATE.} For completeness, we briefly restate the hierarchical learning strategy from Section (ref), now adopted to construct a valid posterior $\Pi_\Delta$ for $\Delta$. Building on the motivation and derivations in Section (ref), construction proceeds by first drawing $b \sim \Pi_b$ using the joint posterior factorization in (ref). Conditional on $b$, the posterior $\Pi_{\Delta \mid b}$ is obtained through the conditional likelihood formulation with a suitably chosen prior; see Section (ref) for details. The {\it marginal posterior for $\Delta$} is then defined as:

where the second equality follows from (ref). This procedure then yields a {\it valid} posterior for $\Delta$, integrating the joint bias information and the hierarchical learning framework, both of which are central to the generalized DRDB procedure.

\paragraph{Posterior aggregation.} Finally, we construct the final posterior $\Pi_\Delta^\mathrm{CF}$ for $\Delta$ using the entire data $\mathcal{D}$, to recover the efficiency lost. Following the DRDB with CF procedure in Section (ref), we combine the posteriors $\Pi_\Delta^{(1)}, \cdots \Pi_\Delta^{(\mathbb{K})}$ obtained from the corresponding splits $(\mathcal{D}_k, \mathcal{D}_k^\-)$ via a CMC approach. We draw independent samples $\{\Delta^{(k)}\}_{k = 1}^\mathbb{K}$ from these posteriors, and define a new random variable:

The {\it aggregated posterior} $\Pi_\Delta^\mathrm{CF}$, {\it our final output}, integrates information on $\Delta$ across all splits, providing a principled and computationally efficient basis for final inference on the ATE (see Theorem (ref)). The main steps of the DRDB procedure with CF are summarized in Algorithm (ref).

Theory

This section develops the theoretical foundations of the DRDB procedure and establishes posterior consistency and BvM–type results (Theorems (ref)--(ref)) for the final DRDB posteriors, under mild regularity conditions on the nuisance parameters. We first present the posterior construction details for DRDB in Section (ref), followed by the result for $\mu_1^\dagger$ in Theorem (ref), a problem of independent interest in missing data theory, and thereafter the main result for the ATE in Theorem (ref).

Likelihood constructions and posterior calculations

We provide a general characterization of the likelihood construction and prior specification used to obtain the posterior of the bias $b(\underbar{m}_t, \underbar{r}_t)$ for $t = 0,1$ and the conditional posterior calculation used in deriving the marginal posterior for $\Delta$ (and $\mu_1$) as discussed in Sections (ref) and (ref).

To avoid repetition, we present a unified procedure applicable to any bias term defined in Sections (ref) and (ref). Likewise, the conditional posterior derivation is also framed generally, covering the computation of conditional posteriors for $\Delta$ and $\mu_1$, given the corresponding bias(es). This construction includes the specific forms used in Sections (ref) and (ref) as special cases.

\paragraph{Bias modeling via targeted modeling strategy.} For notational convenience, we use $\mathcal{N}(\mu, \sigma^2)$ for a Normal distribution with mean $\mu$ and variance $\sigma^2$, and $t_\nu(\eta, c^2)$ for a $t$-distribution with degrees of freedom $\nu >0$, center $\eta$ and scale $c$. Define $b_t := b(\underbar{m}_t, \underbar{r}_t)$, $W\!(\mathbf{Z},\underbar{r}_t,\underbar{m}_t):= \underbar{r}_t(\mathbf{X})\{Y -\underbar{m}_t(\mathbf{X})\}$ and $\sigma^2_t := \mathrm{Var}_{(Y,\mathbf{X}) \in S_t}\{W\!(\mathbf{Z},\underbar{r}_t,\underbar{m}_t)|\underbar{r}_t, \underbar{m}_t\}$ for $t = 0,1$ and $\mathbf{Z} = (Y, \mathbf{X}) \in S (\perp \!\!\! \perp (\underbar{\it m}_t, \underbar{\it r}_t))$. Then, given $\underbar{m}_t \sim \Pi_{m_t}$ and $\underbar{r}_t \sim \Pi_{r_t}$, for $\mathbf{Z}_i \in S_t \subset S$, the weighted observables $W\!(\mathbf{Z}_i,\underbar{r}_t, \underbar{m}_t)$ are i.i.d. with mean $b_t$ and variance $\sigma^2_t$. This motivates a natural {\it working} model based on a Normal distribution with unknown variance. For simplicity, we recommend using an improper prior on the model parameters, though more general priors yield the same asymptotic properties. Let $\mathcal{I}_t$ denote the index set of $S_t$. The model and prior formulation are then given as: for each $t \in \{0,1\}$,

\paragraph{Conditional posterior construction.} To generalize the conditional posterior derivation, we introduce the generic random variables $\theta$ and $\lambda$, where $\theta$ represents either $\Delta$ or $\mu_1$, and $\lambda$ denotes the corresponding bias, i.e., $b_1$ in Section (ref) or $b$ in Section (ref). Since the likelihood is constructed using the entire data $S$ and does not depend on specific properties of the bias term or the target parameter, we are justified in adopting this unified notation.

Let $\varphi(\cdot)$ denote a generic regression function, corresponding to $m_1(\cdot)$ in Section (ref) or $m(\cdot)$ in Section (ref). Given posterior samples $\underline{\varphi} \sim \Pi_\varphi \equiv \Pi_\varphi(\cdot; S^\-)$ and $\underline{\lambda} \sim \Pi_{\lambda} \equiv \Pi_{\lambda}(\cdot; S)$, we have the i.i.d. replicates $\{\underline{\varphi}(\mathbf{X}_i)\}_{i \in \mathcal{I}} \in S$ {\it targeting} $\theta - \underline{\lambda}$ through their mean. Within the target-specific modeling strategy, $\theta - \underline{\lambda}$ can be viewed as a functional of the distribution of $S$, characterized by the summary statistic (mean) of $\underline{\varphi}(\mathbf{X})$. Given the independence of these observables, it is natural to adopt a Normal {\it working} model with unknown variance and place an improper prior (for simplicity again, though more general priors are allowed) on its parameters, yielding an analytically tractable posterior. Specifically, let $\sigma^2_2 := \mathrm{Var}_{\mathbf{X} \in S}\{\underline{\varphi}(\mathbf{X}) | \underline{\varphi}\}$. Then, the resulting model is specified as:

In particular, setting $\theta = \mu_1, \underline{\lambda} = \underbar{b}_1$ and $\underline{\varphi} = \underbar{m}_1 \sim \Pi_{m_1}$ gives the conditional posterior $\Pi_{\mu_1|b_1}$ for $\mu_1$ given $\underbar{b}_1 \sim \Pi_{b_1}$, as in Section (ref). Similarly, setting $\theta = \Delta, \underline{\lambda} = \underbar{\it b}$ and $\underline{\varphi} = \underbar{\it m} \sim \Pi_{m}$ yields the conditional posterior $\Pi_{\Delta|b}$ for $\Delta$ given $\underbar{\it b} \sim \Pi_{b}$, as in Section (ref).

Main results

This section presents the theoretical properties of the proposed DRDB procedure, providing results for both the ATE $\Delta^\dagger$ and the one-arm parameter $\mu_1^\dagger\equiv\mu^\dagger(1) = \mathbb{E}[Y(1)]$.

Let $P$ and $Q$ be two probability measures on a measurable space $(\Omega, \mathcal{B})$. Then, the total variation distance between $P$ and $Q$ is defined as: $d_\mathrm{TV}(P, Q):= \sup_{B \in \mathcal{B}}|P(B) - Q(B)|$.

A direct consequence of Theorems (ref) and (ref) is that DRDB provides natural Bayesian point estimators for $\mu^\dagger(1)$ and the ATE $\Delta^\dagger$ through the {\it posterior means:} $\widehat{\mu}_1(\underbar{m}_1, \underbar{r}_1)$ and $\widehat{\Delta}(\underbar{m}, \underbar{r})$, respectively. In Corollary (ref), we rigorously characterize the theoretical properties of these DRDB point estimators.

Similar asymptotic properties to those in Remark (ref) hold for $\widehat{\mu}_1(\underbar{m}_1, \underbar{r}_1)$, relevant for mean estimation of missing outcomes under MAR tsiatis2007semiparametric. The details are analogous and omitted for brevity.

Comparison with alternative Bayesian debiasing strategies

A main challenge in Bayesian semiparametric inference is that regularization bias from flexible nuisance models can propagate into the posterior for a low-dimensional target parameter, such as the ATE, compromising inferential validity bickel2012semiparametric, castillo2015bernstein. To mitigate this nuisance-induced bias, two prominent strategies have emerged: {\it prior modification}, which tailors prior specification to the semiparametric model structure ray2020semiparametric, breunig2025double; and {\it posterior correction}, which applies a post-hoc adjustment using the efficient influence function (EIF) yiu2025. A complementary approach by luo2023semiparametric constructs posteriors via exponentially tilted empirical likelihood and establishes BvM results for partially linear and parametric models. Our DRDB procedure introduces an alternative perspective by embedding debiasing directly into the modeling process through {\it targeted learning}. Below, we compare DRDB with these two state-of-the-art alternatives, focusing on the seminal works of ray2020semiparametric and yiu2025, both mainly interested in the mean outcome under MAR, closely related to the one-arm case for our ATE setting, discussed in Section (ref).

The prior modification approach proposed by ray2020semiparametric models the {\it full data distribution} with nonparametric priors, innovatively augmenting the prior for the outcome regression with an estimator of the propensity score (PS) obtained from an independent auxiliary data. This augmentation perturbs the prior in the model’s least favorable direction of the semiparametric model, and thereby mitigates nuisance-induced bias. While theoretically elegant, this approach requires customized prior design and strong conditions, including Donsker assumptions and smoothness constraints, limiting the use of machine learning or high dimensional nuisance models. Moreover, its theoretical guarantees require both nuisance models to be correctly specified.

The posterior correction approach in yiu2025 takes a different route by adding a stochastic correction to posterior draws, inspired by the frequentist {\it one-step estimator} van2000asymptotic. Using the EIF and the {\it Bayesian bootstrap} rubin1981bayesian, it projects draws towards the truth along the most informative direction, achieving bias reduction. While this approach can target multiple functionals, it relies on posterior draws for the full data model, and, similar to DRDB, its theoretical guarantees rely on product-type rate conditions on the nuisance posteriors. In addition, it requires Donsker-type conditions van2000asymptotic, which are often restrictive for flexible or high-dimensional nuisance models, whereas DRDB avoids the explicit need for such conditions (see Assumption (ref)) via its distinct use of cross-fitting.

Cross-fitting (CF) is a well-established tool in the frequentist literature, commonly used to relax Donsker-type conditions on nuisance parameters chernozhukov2018double. In DRDB, however, CF is {\it not} merely a technical device; it is a crucial component of our debiasing mechanism. Within the DRDB framework, CF is methodologically essential (see Remark (ref)), and also underpins a novel {\it Bayesian analogue}, an aggregation strategy that combines {\it posteriors} across splits (folds), enabling principled Bayesian semiparametric inference while leveraging the full data efficiently.

DRDB departs from both methods in philosophy and implementation. Rather than modifying the prior or correcting the posterior, DRDB embeds debiasing within the modeling process through the {\it debiased representation}. It explicitly identifies and {\it learns the bias as a separate target}, using summary statistics that are directly informative about the ATE and its bias. Furthermore, the role of the propensity score underlines these differences: it enters externally to guide prior design, while in DRDB it arises {\it naturally} through density ratio weighting. Likewise, the use of independent data differs: ray2020semiparametric primarily leverages it for technical convenience, whereas DRDB uses it to {\it validate} the debiased representation and strengthen bias correction.

Theoretically, DRDB requires only {\it high-level posterior contraction} for the nuisance models, and mild moment conditions, making it compatible with flexible nuisance models. Most importantly, DRDB achieves {\it Bayesian double robustness}: ATE posterior remains consistent and contracts at the rate of the well-specified nuisance, even if the other is misspecified. This guarantee is stronger than those of both prior augmentation and one-step posterior correction methods, offering greater stability under model misspecification. Computationally, DRDB is also far simpler and more scalable. Whereas yiu2025 require a full set of posterior draws for all nuisance parameters, DRDB needs only a {\it single} posterior draw per nuisance per cross-fitting fold. Compared with ray2020semiparametric, which demands intricate prior customization, DRDB’s modular structure allows seamless use of standard Bayesian regression tools without model-specific tuning.

A further distinction of DRDB is the generality of its methodological framework, which extends seamlessly beyond the ATE to a broad class of causal estimands. By representing each estimand as a weighted functional and adapting the debiasing and retargeting steps accordingly, DRDB maintains inferential validity and computational scalability under these broader settings. The detailed formulation of this extension is presented in Section (ref) of the \hyperref[sec_supplementary]{Supplementary Material}.

DRDB unifies the theoretical strengths of prior modification and posterior correction while also introducing new modeling perspectives and methodological advances that extend its applicability. By embedding debiasing directly into the modeling process, it enables flexible nuisance estimation, ensures valid inference under mild conditions, and remains computationally efficient. These features establish DRDB as a robust, theoretically grounded, and practically scalable framework for Bayesian causal inference and, more broadly, Bayesian semiparametric inference in general.

Numerical studies

We evaluate the finite-sample performance of the proposed DRDB procedure for both estimation and inference of the ATE through extensive simulation studies across various data-generating mechanisms and Bayesian methods for nuisance estimation, including both well-specified and misspecified settings. The mean of the DRDB posterior $\Pi_\Delta^\mathrm{CF}$ serves as our point estimator. We report the empirical bias ({\bf Bias}) and mean squared error ({\bf MSE}) for estimation accuracy. For inference evaluation, we report the empirical coverage probabilities ({\bf Cov}) and average lengths of the 95% credible intervals ({\bf CI-Len}) based on 1000 posterior samples from $\Pi_\Delta^\mathrm{CF}$. For the number of folds, we set $\mathbb{K} = 5$ for computational efficiency. All reported results are based on 500 replications. We study two scenarios: one where {\it both} nuisance models are well-specified (Section (ref)), and another where {\it only one} of them is well-specified (Section (ref) of the \hyperref[sec_supplementary]{Supplementary Material}).

The following notations are used throughout this section. For any integer $p \geq 1$ and $v \in \mathbb{R}$, let $v_p$ be the vector $v_p := (v, \dots, v)' \in \mathbb{R}^{p \times 1}$. Let $I_p$ denote the $p \times p$ identity matrix, $\mathcal{N}_p(\mu_p, \Sigma_p)$ denotes the $p$-variate Gaussian distribution with mean vector $\mu_p \in \mathbb{R}^{p}$ and covariance matrix $\Sigma_p \in \mathbb{R}^{p \times p}$.

Simulation results

Throughout, we set $n = 1000$ and consider $p = 10$, $50$, and $200$, representing low, moderate, and high-dimensional settings, respectively. For each $i = 1, \dots, n$, the covariate vector is generated as: $\mathbf{X}_i \overset{\text{i.i.d.}}\sim \mathcal{N}_p(0_p, I_p)$. Conditional on $\mathbf{X}_i$, the treatment assignment follows: $T_i| X_i \sim \mathrm{Ber}\{e^\dagger(\mathbf{X}_i)\}$, where $e^\dagger(\mathbf{X}_i) = 1/\{1 + \exp^{-(\mathbf{X}_i'\beta_3 - 0.08)}\}$ and $\beta_3 = (0.35_2, 0_{p-2})$, ensuring the positivity condition in Assumption (ref). Given $\mathbf{X}_i$, the potential outcomes are generated as: $Y_i(t) \sim \mathcal{N}(m_t^\dagger(\mathbf{X}_i), \sigma_t^2)$, with $m_1^\dagger(\mathbf{X}) = 5 + 2 \mathbf{X}'\beta_1$ and $m_0^\dagger(\mathbf{X}) = 3 + \mathbf{X}'\beta_0$, and variances $\sigma_t^2 = \mathrm{Var}\{m_t^\dagger(\mathbf{X})\}/5$ for $t = 0, 1$. The observed outcome is therefore $Y_i = T_i Y_i(1) + (1-T_i) Y_i(0)$, and the observed data is $\mathcal{D} = \{(Y_i, \mathbf{X}_i, T_i)\}_{i=1}^n$. The regression coefficients $\beta_1 = \beta_0$ are set as $ (1_{s/2}, 0.5_{s/2}, 0_{p-s})$, where $s$ denotes {\it sparsity}. For $p = 10$, we use $s = 3$ and $s = 10$; for $p = 50$ and $p = 200$, we take $s \approx \sqrt{p}$ and $s \approx p/4$ to represent sparse and moderately dense regimes.

For illustrative purposes, we employ the sample mean $\widehat{p}_1: = n^{-1}\sum_{i = 1}^n T_i$ as a point estimator for $p_1 := \mathbb{P}(T = 1)$. For the posterior $\Pi_e$, we only use sparse Bayesian ({\tt BS}) logistic regression with nonlocal priors (NLP) johnson2012bayesian for simplicity. For the posteriors $\Pi_{m_1}$ and $\Pi_{m_0}$, in addition to {\tt BS} linear regression with NLP, we consider Bayesian ridge regression ({\tt BR}) and {\tt BART} bart2010, implemented using the R package BART. For parametric methods, we consider the Gaussian linear and logistic regression working models: for $i = 1, \dots, n$, $Y_i | \mathbf{X}_i, T_i = 0, \gamma_0, \theta_0,\sigma \overset{\text{i.i.d.}}\sim \mathcal{N}(\gamma_0 + \mathbf{X}_i'\theta_0, \sigma^2)$, $Y_i| \mathbf{X}_i, T_i = 1, \gamma_1, \theta_1,\tau \overset{\text{i.i.d.}}\sim \mathcal{N}(\gamma_1 + \mathbf{X}_i'\theta_1, \tau^2)$ and $T_i | \mathbf{X}_i \sim \mathrm{Ber}\{e(\mathbf{X}_i)\}$, where $e(\mathbf{X}_i) = 1/\{1+ \exp^{-(\gamma_3 + \mathbf{X}_i'\theta_3)}\}$. For {\tt BR}, we employ a Gaussian prior on the regression coefficients and an improper prior on the variance parameter. The ridge parameter $\lambda$ is estimated using an empirical Bayes approach, with the point estimate $\widehat{\lambda}$ obtained via the R package glmnet. For {\tt BS}, posterior samples for $(\gamma_0, \theta_0)$, $(\gamma_1, \theta_1)$ and $(\gamma_3, \theta_3)$ are obtained using the R package mombf. Finally, as a performance benchmark, we report the results for the frequentist {\it oracle} estimator, constructed using the empirical mean of the EIF of the ATE hahn98 over $\mathcal{D}_n$ using the {\it true} nuisance parameters, denoted as {\tt Oracle}.

Table (ref) reports the estimation and inference results for the ATE. Across all scenarios and nuisance estimation methods, DRDB performs nearly identically to the {\tt Oracle} estimator. In both low- and moderate-dimensional settings ($p = 10, 50$), and across sparse ($s = \sqrt{p}$) and moderately dense ($s = p/4$) regimes, all DRDB variants yield negligible differences from {\tt Oracle}. Even DRDB-B performs comparably well, though with slightly higher finite-sample bias and wider credible intervals, reflecting the slower convergence of nonparametric nuisance estimation. As the dimensionality increases to $p = 200$, the role of sparsity becomes more pronounced. The sparsity-adaptive DRDB-S {\it continues} to perform close to {\tt Oracle}, effectively leveraging the underlying sparse structure, whereas DRDB-R and DRDB-S exhibit increased finite-sample bias and wider CIs, with DRDB-B most affected due to slower convergence inherent to nonparametric procedures.

Table (ref) shows that DRDB {\it consistently} achieves coverage probabilities near the {\it nominal} 95% level across nearly all settings with various dimensions, sparsity levels, and the nuisance estimation methods. In high-dimensional scenarios, DRDB-B tends to produce slightly conservative coverage and wider CIs, reflecting the effect of slower convergence rates and increased finite-sample bias associated with higher parameter dimensionality. Nevertheless, CI lengths {\it remain comparable} to those of the {\tt Oracle}, and DRDB-S consistently {\it maintains} coverage near the nominal level. These patterns highlight DRDB’s {\it robustness} and {\it adaptability}: when nuisance models capture the true underlying structure, its performance remains near-optimal {\it even} in high-dimensional regimes.

Figures (ref) and (ref) visually support these findings. Across all scenarios, DRDB posteriors exhibit approximately {\it Gaussian shapes}, always concentrated around the truth $\Delta^\dagger = 2$, with posterior means (from box plots) {\it tightly} centered near the truth. In low or moderate-dimensional regimes ($p = 10$ or $p = 50$), all DRDB variants produce nearly {\it identical} posteriors that closely {\it align with} the Oracle, confirming the {\it stability} and {\it efficiency} of DRDB when nuisance parameters are accurately estimated. In the high-dimensional settings ($p = 200$), however, the differences between nuisance models become more pronounced: DRDB-S maintains the most concentrated posterior, nearly matching the Oracle, while DRDB-R and DRDB-B show wider and heavier-tailed posteriors, reflecting slower nuisance convergence and the inherent difficulty of ridge and nonparametric estimators in such settings. Yet, all density curves {\it remain} well centered around $\Delta^\dagger = 2$, and the credible interval lengths are largely comparable to the Oracle’s. Overall, the plots confirm DRDB’s robust finite-sample performance, while clearly revealing how high dimensionality amplifies the effect of nuisance model choice on second-order properties of posterior concentration.

Finally, while our analysis herein focuses on correctly specified models, high-dimensional settings inherently introduce a {\it soft} form of misspecification due to finite-sample nuisance estimation bias. Even in such cases, DRDB {\it consistently} estimates the ATE across nuisance methods, demonstrating its double robustness: as long as the product of nuisance convergence rates exceeds the parametric rate, DRDB yields near-Oracle estimates with {\it accurate coverage}. Under explicit functional-form misspecification (see Section (ref) in the \hyperref[sec_supplementary]{Supplementary Material}), DRDB continues to produce stable estimates across all nuisance models and maintains {\it valid coverage} with only mildly wider credible intervals. These findings support its {\it double robustness}, and highlight a distinct {\it advantage} of the {\it Bayesian} framework as well. Overall, the results confirm that DRDB offers robust, efficient, and {\it reliable inference} across diverse settings, highlighting its {\it insensitivity} to nuisance estimation.

Real data application

We apply the proposed DRDB approach to evaluate {\it the effect of smoking cessation on weight gain} using data from the National Health and Nutrition Examination Survey Epidemiologic Follow-up Study (NHEFS). The NHEFS is a longitudinal study initiated by the National Center for Health Statistics and the National Institute on Aging, in collaboration with other agencies of the U.S. Public Health Service Hernan2020, and is well-studied in the causal inference literature Ertefaie2022, zhang2023double. Our analysis focuses on a subset of $n = 1566$ individuals who were cigarette smokers aged 25--74 years at baseline in 1971, and had a follow-up visit in 1982. The treatment $T \in \{0,1\}$ ({\tt qsmk}) equals 1 if the individual quit smoking before the follow-up and 0 otherwise. The outcome $Y \in \mathbb{R}$ ({\tt wt82_71}) is the weight gain (in kg), calculated as the difference between body weight at the follow-up and baseline weight. We use the same nine covariates as in Ertefaie2022: {\tt sex}, {\tt race}, {\tt age}, {\tt education}, {\tt smokeintensity}, {\tt smokeyrs}, {\tt exercise}, {\tt active}, and {\tt wt71}. A detailed data description is available at \url{https://miguelhernan.org/whatifbook}.

Our goal is to estimate the ATE of the treatment on body weight gain. As a baseline for comparison, we include the naive estimator ({\tt Naive}), defined as the difference in empirical mean outcomes between the treatment groups tsiatis2007semiparametric. Using the Bayesian nuisance estimation methods given in Section (ref), we compute the DRDB posterior of $\Delta$ via Algorithm (ref) with $\mathbb{K} = 5$. For each posterior, we draw 1000 samples to compute Monte Carlo approximations of the posterior mean (as the point estimate) and the 2.5% and 97.5% quantiles to construct 95% credible intervals (CIs).

Table (ref) reports the ATE estimates, the 95% CIs, and the respective CI lengths (CI-Leng). All methods yield a positive ATE, indicating that quitting smoking is associated with weight gain. A notable disparity exists between the naive estimate of 2.541 and the DRDB-based estimates, which are consistently higher and clustered in the 3.31 to 3.54 range. This gap strongly suggests the presence of confounding (given the observational setting), which the naive estimator fails to address. Further, our estimates {\it align} closely with those of Ertefaie2022, who reported ATE estimates between 3.20 and 3.42 using IPW estimators on the same data, providing strong external validation for our method. Moreover, DRDB yields substantially narrower CIs, with lengths of approximately 1.55, compared with the naive estimator (length $= 1.91$) and the IPW-based CIs in Ertefaie2022 (length $\approx 2.34$ -- $2.41$), demonstrating {\it improved efficiency}. This consistency in results across various nuisance estimation methods highlights the {\it robustness} of DRDB to the choice of nuisance model. Overall, these findings indicate the presence of notable confounding via $\mathbf{X}$, and support a causal effect of smoking cessation on weight gain after adjusting for confounders.

Concluding discussion

We proposed a DRDB procedure for estimating {\it causal functionals} within a Bayesian framework, focusing on the ATE and the single-arm mean $\mu^\dagger(1)$. DRDB builds on {\it two key ideas}: (i) explicit separation of the nuisance estimation from inference on the target, through a {\it Bayesian debiasing} mechanism (coupled with {\it data splitting}); and (ii) a {\it targeted modeling} strategy, along with {\it hierarchical learning} and use of {\it posterior aggregation}. The debiasing step corrects nuisance-induced bias via a density-ratio–based retargeting, while targeted modeling ensures efficient use of relevant summaries. All these aspects put together make DRDB both theoretically robust/efficient and computationally scalable. Theorems (ref)--(ref) establish BvM theorems for DRDB with matching frequentist guarantees when both nuisances are correctly specified, as well as a Bayesian analogue of the frequentist double robustness, while operating within the posterior framework. Computationally, DRDB is flexible, allowing for {\it off-the-shelf} nuisance model choices, admits a simple final posterior, and requires only a {\it single} posterior draw from each nuisance per CF fold, ensuring its scalability in high dimensions. Finally, our framework also readily extends to a wide class of policy-relevant causal estimands, as shown in Section (ref) of the \hyperref[sec_supplementary]{Supplementary Material}, with similar theoretical guarantees expected to hold therein as well. Overall, DRDB provides a {\it novel} Bayesian perspective on debiasing and double robustness, offering both theoretical guarantees and practical feasibility, and opens the door to broader applications of Bayesian targeted semiparametric learning.

\phantomsection \addcontentsline{toc}{section}{Acknowledgments and funding}

Acknowledgments and funding

The authors gratefully acknowledge funding support from the NSF grants: NSF-DMS 2113768 (Abhishek Chakrabortty) and NSF-DMS 2210689 (Anirban Bhattacharya) towards partial support of this research.

\phantomsection \addcontentsline{toc}{section}{Data availability}

Data availability

The dataset used in this paper is publicly available. The code used to generate the simulation results is available at \url{https://github.com/gozdesert/DRDB}.

\phantomsection \addcontentsline{toc}{section}{Supplementary material}