# Semiparametric posterior corrections

Andrew Yiu, Edwin Fong, Chris Holmes, Judith J. Rousseau

► To cite this version:

Andrew Yiu, Edwin Fong, Chris Holmes, Judith J. Rousseau. Semiparametric posterior corrections. Journal of the Royal Statistical Society: Series B, 2025, 87 (4), pp.1025-1054. 10.1093/jrsssb/qkaf005. hal-05354491

HAL Id: hal-05354491

https://hal.science/hal-05354491v1

Submitted on 7 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Semiparametric posterior corrections

Andrew Yiu<fnref n="1" />, Edwin Fong<fnref n="2" />, Chris Holmes<fnref n="1" />,<fnref n="3" /> and Judith Rousseau<fnref n="1" />,<fnref n="4" />

*Address for correspondence:* Andrew Yiu, School of Mathematical Sciences, Building 54, Highfield Campus, Southampton SO14 3ZH, UK. Email: a.h.m.yiu@soton.ac.uk

## Abstract

We present a new approach to semiparametric inference using corrected posterior distributions. The method allows us to leverage the adaptivity, regularization, and predictive power of nonparametric Bayesian procedures to estimate low-dimensional functionals of interest without being restricted by the holistic Bayesian formalism. Starting from a conventional posterior on the whole data-generating distribution, we correct the marginal posterior for each functional of interest with the help of the Bayesian bootstrap. We provide conditions for the resulting *one-step posterior* to possess calibrated frequentist properties and specialize the results for several canonical examples: the integrated squared density, the mean of a missing-at-random outcome, and the average causal treatment effect on the treated. The procedure is computationally attractive, requiring only a simple, efficient postprocessing step that can be attached onto any arbitrary posterior sampling algorithm. Using the ACIC 2016 causal data analysis competition, we illustrate that our approach can outperform the existing state-of-the-art through the propagation of Bayesian uncertainty.

**Keywords:** Bayesian inference, causal inference, double robustness, semiparametric theory

## 1 Introduction

Functional estimation is emerging as a major paradigm in modern statistics (Buja et al. 2019a, 2019b; van der Laan & Rose, 2011; Vansteelandt & Dukes, 2022). The definition of the estimand—the target of inference—is placed at the forefront of the statistical analysis, requiring the user to carefully articulate their scientific objectives. This helps to avoid the common pitfall in which the choice of estimand is informed by the choice of statistical model, rather than the other way round. Causal inference, for example, has a natural affinity with functional estimation through the nonparametric identification of causal effects (Hernán & Robins, 2020; Pearl, 2009).

An important advantage of the functional approach is that we are free to estimate low-dimensional functionals without being restricted to low-dimensional models. Nonparametric methods provide greater flexibility in capturing complex relationships between variables, mitigating the risks of model misspecification bias. In particular, it is becoming increasingly popular to leverage highly flexible black-box algorithms (e.g. random forests, gradient boosting, deep learning) to estimate nuisance parameters (Chernozhukov et al., 2018; van der Laan & Rubin, 2006; Vansteelandt & Dukes, 2022).

Using infinite-dimensional models to estimate low-dimensional estimands is often referred to as *semiparametric inference*. By exploiting the smoothness (or *differentiability*) of the target functional, it may be possible to obtain estimators with parametric properties like $\sqrt{n}$-consistency and asymptotic normality, even if the nuisance parameter estimators converge at a slower rate (Hines et al., 2022; Murphy & van der Vaart, 2000; van der Laan & Rose, 2011). Some authors

<footnote n="1">¹Department of Statistics, University of Oxford, Oxford, UK</footnote>

<footnote n="2">²Department of Statistics and Actuarial Science, University of Hong Kong, Hong Kong</footnote>

<footnote n="3">³Nuffield Department of Medicine, University of Oxford, Oxford, UK</footnote>

<footnote n="4">⁴Ceremade, CMRS, UMR 7534, Université Paris-Dauphine, PSL University, Paris, France</footnote>
have suggested that semiparametric methods can bridge the ‘two cultures’ of statistical modelling (Breiman, 2001) by combining the interpretability and efficiency of parametric inference with the robustness of black-box algorithms (Kennedy et al., 2021; Ogburn & Shpitser, 2021; Vansteelandt, 2021).

Nonparametric Bayesian approaches such as Bayesian additive regression trees (BART) (Chipman et al., 2007) have performed especially well in empirical studies for both prediction and inferring low-dimensional estimands such as causal effects (Dorie et al., 2019; Hahn et al., 2020; Hill, 2011). For parametric Bayesian models, the influence of the prior disappears asymptotically to first-order thanks to the renowned Bernstein-von Mises (BvM) theorem (Le Cam, 1986; van der Vaart, 1998). For infinite-dimensional models, however, the impact of the prior may remain significant in the asymptotic regime, and the prior specification cannot be completely justified by subjective beliefs (Diaconis & Freedman, 1986; Rousseau, 2016). This has led to a rich literature on the frequentist properties of Bayesian nonparametric methods, with a strong emphasis on posterior contraction rates (e.g. S. Ghosal & van der Vaart, 2017). In particular, many Bayesian approaches have been shown to enjoy excellent properties in terms of adaptive posterior contraction rates due to the flexibility enabled by hierarchical prior modelling (e.g. S. Ghosal & van der Vaart, 2017; Rousseau, 2016).

Bayesian approaches naturally deal with nuisance parameters by integrating them out, which allows for coherent inference on possibly multiple parameters of interest, whether they are estimated simultaneously or sequentially (e.g. J. O. Berger et al., 1999). Unfortunately, good adaptive contraction rates do not necessarily translate to good behaviour for marginal posterior distributions of specific functionals, even to first-order. Although there has been a growing literature on the existence of BvM theorems in the semiparametric setting (Castillo, 2012; Castillo & Nickl, 2014; Castillo & Rousseau, 2015; Rivoirard & Rousseau, 2012), the results have been only partially positive, and it is now recognized that a given prior will perform well for some functionals of interest but not for others. This is perhaps unsurprising given that Bayesian inference is a ‘plug-in’ approach, and similar phenomena have been observed for semiparametric plug-in estimators in the frequentist literature (Bickel & Ritov, 2003; Robins et al., 2017; van der Laan & Rubin, 2006). The reason is that flexible nonparametric methods must employ regularization to achieve good performance, and the ‘regularization bias’ can bleed into the corresponding plug-in estimator, precluding fast convergence rates (Chernozhukov et al., 2018; van der Laan & Rubin, 2006; van der Vaart, 2014).

For certain examples, it has been demonstrated that the asymptotic posterior bias for a non-parametric Bayesian model can be controlled by tailoring the prior specification to the estimand (Castillo & Rousseau, 2015; Ray & van der Vaart, 2020). This is undesirable from a practical perspective, however, as it will require careful tuning (or even rewriting) of any existing software for posterior computation. Moreover, a different modification will be required for each estimand of interest. In any case, studying the existence of a BvM theorem remains a formidable open question for many popular families of priors, such as nonparametric mixtures. For the implementations of nonparametric Bayesian procedures that are widely applied in practice, existing empirical evidence suggests that the lack of a semiparametric BvM theorem represents the rule, rather than the exception (Dorie et al., 2019; Hahn et al., 2020; Ray & Szabó, 2019).

To address these issues, we introduce a simple postprocessing procedure that starts from a given posterior distribution on the whole data-generating distribution and then corrects the marginal posterior for each functional of interest. This operates by adding a stochastic term based on the *efficient influence function* of the functional (Pfanzagl, 1982; van der Vaart, 1991, 1998) and the Bayesian bootstrap (Rubin, 1981), which plays a crucial role for correcting not only the bias of the posterior but also its shape. The original user-specified Bayesian model can be left untouched, acting as a central hub from which we can target each estimand of interest individually.

We provide general conditions for these new *one-step posteriors* to satisfy a BvM theorem, which ensures that central credible regions achieve approximately nominal coverage with the semiparametric efficient size. Unlike the semiparametric BvM theorems derived for standard nonparametric posteriors (e.g. Castillo, 2012; Castillo & Rousseau, 2015), our conditions do not depend delicately on the likelihood and prior. Instead, we require only contraction rates and complexity bounds on the posterior with empirical process theory, leading to a Bayesian counterpart to the classical theory of influence-function-based estimation (Hines et al., 2022;
Pfanzagl, 1982; van der Laan & Robins, 2003; van der Laan & Rose, 2011). In particular, the one-step posterior could be interpreted as a natural analogue of the *one-step estimator* (Newey et al., 1998; Pfanzagl, 1982) for correcting an entire posterior distribution rather than just a point estimator. Moreover, we can use the one-step posterior to provide insight on the potential presence of regularization bias in the initial posterior, which opens up an interesting direction for performing sensitivity analysis.

Our methodology, which we introduce in Section 2, is computationally efficient and attaches onto any existing posterior sampling implementation without modification of the original algorithm. The procedure takes each posterior sample and adds a randomized correction term, which only requires drawing an independent set of uniform Dirichlet weights of lengths equal to the sample size. A single set of posterior samples can be retained for simultaneous inference of multiple functionals of interest.

We apply our approach to the classic example of estimating the integrated squared density in Section 3. This was previously studied in a Bayesian setting by Castillo and Rousseau (2015) but only for random histogram priors on the interval [0, 1]. We verify the conditions for the more complex class of Dirichlet process Gaussian location mixtures (Shen et al., 2013). To the best of our knowledge, this is the first BvM theorem associated with priors based on nonparametric mixture models.

In Section 4, we study the estimation of the mean of an outcome that is missing-at-random. This is a problem that has received much attention due to its connections with estimating the average treatment effect in causal inference (Kang & Schafer, 2007; Robins et al., 2017). The conditions for the BvM theorem are expressed here in terms of the propensity score and outcome regression models, showing that the asymptotic behaviour of the one-step posterior depends crucially on the combined contraction rates of both parameters. For this specific example in the binary outcome setting, Ray and van der Vaart (2020) proposed a correction approach that has similar motivations and asymptotic considerations to ours, but it requires a modification to the prior that can be relatively complicated to implement. Working within their set-up, we study the behaviour of the one-step posterior and provide a comparison of theoretical assumptions.

The last example is the average treatment effect on the treated (Heckman & Robb, 1985; Rubin, 1977), which has not been studied before from a large-sample Bayesian perspective. Due to its complex structure, this estimand poses some theoretical challenges. In Section 5, we introduce a slight modification of the correction algorithm, taking inspiration from estimating equation methodology. Extensions of this approach to sample conditional treatment effects can be found in online supplementary material, Section S6.

We evaluate our methodology empirically in Section 6. This includes a large-scale comparison with state-of-the-art causal algorithms using the ACIC 2016 data analyses competition (Dorie et al., 2019). Our method—applied to BART (Chipman et al., 2007)—is the best performer in terms of both point and interval estimation. In fact, no other method was able to obtain nominal coverage across all of the data-generating mechanisms. We suggest that this superior performance could be attributed to the propagation of Bayesian uncertainty through the posterior correction to regularize the estimation, which cannot be exploited by the classical semiparametric approaches. Thus, we argue that the one-step posterior correction should appeal not only to Bayesians, but also to classically minded users interested in leveraging the power of Bayesian algorithms while obtaining frequentist-calibrated inference.

# 2 Methodology

## 2.1 Set-up and background

Suppose that we observe independent and identically distributed data $Z_{1:n} = (Z_1, \dots, Z_n)$ from a distribution $P_0$ known to belong to a set of probability measures $\mathcal{P}$ on a Polish sample space $(\mathcal{Z}, \mathcal{A})$. We will use the shorthand notation $P[f] = \int f(z) dP(z)$ for $P \in \mathcal{P}$. In particular, $\mathbb{P}_n[f] = n^{-1} \sum_{i=1}^n f(Z_i)$, where $\mathbb{P}_n$ is the empirical measure, and we let $\mathbb{G}_n = \sqrt{n}(\mathbb{P}_n - P_0)$ denote the empirical process. For $P \in \mathcal{P}$, let $L_2(P)$ be the Hilbert space consisting of all measurable functions $h: \mathcal{Z} \to \mathbb{R}$ with $P[h^2] < \infty$ equipped with the inner product $\langle h_1, h_2 \rangle_P = P[h_1 h_2]$ and norm $\|h\|_P = \sqrt{P[h^2]}$.
In a Bayesian analysis, we equip $\mathcal{P}$ with a $\sigma$-field $\mathcal{B}$ such that $(\mathcal{P}, \mathcal{B})$ is a standard Borel space, and we specify a prior probability distribution $\Pi$ on $P$. Let $\Pi(\cdot | Z_{1:n})$ denote a version of the posterior distribution given $Z_{1:n}$. The target estimand is $\chi(P_0)$, where $\chi: \mathcal{P} \to \mathbb{R}$ is a measurable functional. For notational clarity, we have chosen the range space to be one-dimensional; generalizing our results to higher dimensions is straightforward (e.g. by applying the Cramèr-Wold device: see p. 16 of van der Vaart, 1998). This set-up also includes 'separated' semiparametric problems in which we have $\mathcal{P} = \{P_{\theta, \eta}\}$ and $\chi(P_{\theta, \eta}) = \theta$, where $\theta \in \mathbb{R}^k$ for some $k \ge 1$, and $\eta$ is a possibly infinite-dimensional nuisance parameter.

An important special case is $\mathcal{P} = M(\mathcal{Z})$—the set of all probability measures on $(\mathcal{Z}, \mathcal{A})$—equipped with its Borel $\sigma$-field $\mathcal{B}_M$ for the weak topology (J. Ghosh & Ramamoorthi, 2003). Dirichlet processes (Ferguson, 1973) form the canonical family of prior distributions on $(M(\mathcal{Z}), \mathcal{B}_M)$. If $P \sim DP(\alpha)$ for a finite positive measure $\alpha$ on $(\mathcal{Z}, \mathcal{A})$, then $DP(\alpha + n\mathbb{P}_n)$ is a version of the posterior given $Z_{1:n}$. The process $DP(n\mathbb{P}_n)$ is known as the Bayesian bootstrap (Rubin, 1981), which can be interpreted as the noninformative limiting posterior for any sequence of Dirichlet process models in which the prior information $\alpha$ goes to zero. The Bayesian bootstrap plays a central role in our methodology; we denote it by $\Pi_{BB}(\cdot | Z_{1:n})$.

It is now necessary to introduce some concepts from semiparametric theory. We refer the reader to online supplementary material, Section S1 for a detailed overview of the following definitions. The functional $\chi$ is required to be differentiable at all $P \in \mathcal{P}$ in the sense of van der Vaart (1991). The efficient influence function of $\chi$ at $P$ is a measurable function $\dot{\chi}_P: \mathcal{Z} \to \mathbb{R}$ satisfying $P[\dot{\chi}_P] = 0$ and $\|\dot{\chi}_P\|_P < \infty$. A sequence of regular estimators $\dot{\chi}_n = \dot{\chi}_n(Z_{1:n})$ is said to be asymptotically efficient at $P_0$ if

$$ \sqrt{n}(\dot{\chi}_n - \chi(P_0)) = \sqrt{n}\mathbb{P}_n[\dot{\chi}_{P_0}] + o_{P_0}(1). \quad (1) $$

The resulting limiting distribution $\mathcal{N}(0, \|\dot{\chi}_{P_0}\|_{P_0}^2)$ is optimal in terms of both the convolution theorem (e.g. Theorem 25.20 of van der Vaart, 1998) and the local asymptotic minimax theorem (van der Vaart, 1992).

For smooth parametric models, it is well-known that—under mild regularity conditions—the marginal posterior of the parameter converges to a normal distribution centred at an efficient estimator (e.g. the maximum likelihood estimator) with covariance equal to the inverse Fisher information divided by the sample size. This is the (efficient) BvM theorem (Chapter 10 of van der Vaart, 1998). Consequently, the frequentist coverage of many credible regions—such as balls around the posterior mean, equal tails intervals, highest posterior density regions etc.—will be asymptotically equal to its nominal level, and their sizes will be asymptotically optimal.

A semiparametric BvM theorem describes the marginal posterior of the target functional on an infinite-dimensional model. Similar to the parametric case, we would like the marginal posterior to converge to a normal distribution centred at an efficient estimator with the efficient asymptotic variance. In this article, we study weak convergence of the posterior distributions, which we control using the bounded Lipschitz distance $d_{BL}$ as in Castillo and Rousseau (2015) (see also online supplementary material, Section S2).

**Definition 1** Let $\mathcal{L}_{\Pi}(\sqrt{n}(\chi(P) - \hat{\chi}_n) | Z_{1:n})$ denote the posterior law of $\sqrt{n}(\chi(P) - \hat{\chi}_n)$, where $\hat{\chi}_n$ is the asymptotically efficient sequence $\hat{\chi}_n = \chi(P_0) + \mathbb{P}_n[\dot{\chi}_{P_0}]$. The posterior satisfies the semiparametric BvM theorem if

$$ d_{BL}(\mathcal{L}_{\Pi}(\sqrt{n}(\chi(P) - \hat{\chi}_n) | Z_{1:n}), \mathcal{N}(0, \|\dot{\chi}_{P_0}\|_{P_0}^2)) \xrightarrow{P_0} 0. $$

More informally, we can say that the posterior law of $\sqrt{n}(\chi(P) - \hat{\chi}_n)$ 'converges weakly to $\mathcal{N}(0, \|\dot{\chi}_{P_0}\|_{P_0}^2)$ in probability'.

In the non-Bayesian setting, efficient estimators have been obtained by one-step estimation (Newey et al., 1998; Pfanzagl, 1982), which operates by performing bias corrections that specifically target the estimand. Suppose we have an initial estimate $\hat{P}$ of the data distribution, which induces a plug-in estimate $\chi(\hat{P})$. Since the efficient influence function acts as a first-order
distributional derivative of the functional (Fisher & Kennedy, 2021; Robins et al., 2017), we might expect the ‘second-order remainder’

$$r_2(P_0, P) = \chi(P_0) - \chi(\hat{P}) - P_0[\dot{\chi}_{\hat{P}}]$$

to have some form of quadratic dependence on the estimation error of $\hat{P}$. This suggests that $\chi(\hat{P}) + P_0[\dot{\chi}_{\hat{P}}]$ would be an improvement over the plug-in estimator, but the true data distribution is of course unknown. Replacing $P_0$ with the empirical measure $\mathbb{P}_n$ leads to

$$\hat{\chi}_{1-step} = \chi(\hat{P}) + \mathbb{P}_n[\dot{\chi}_{\hat{P}}].$$

which is called the *one-step estimator*. We review sufficient conditions for the one-step estimator to be asymptotically efficient in online supplementary material, Section S1.2; as we will discuss in Section 2.3, the assumptions for our approach share key conceptual similarities.

**Example 1** (The mean of an outcome that is missing-at-random). Suppose the data takes the form $Z = (X, A, AY)$, where X is a vector of covariates, A is a binary missingness indicator, and Y is the outcome variable of interest that is unobserved when $A = 0$. The target estimand is the outcome mean $\chi(P) = E_P[Y]$, which is identified by

$$\chi(P) = E_P[E_P[Y | A = 1, X]]$$

under the *ignorability* (or *missing-at-random*) (Rubin, 1976) assumption $A \perp Y | X$, and the *positivity* (or *overlap*) assumption: $P(A = 1 | X) > 0$ with $P$-probability 1. In practical terms, the missingness is deemed to be uninformative with respect to the outcome after adjusting for a set of measured covariates. This example is closely related to estimating the average causal effect of a binary treatment in causal inference (see online supplementary material, Section S6 or Morgan & Winship, 2007).

The data distribution $P$ is naturally parameterized in terms of the ‘propensity score’ $\pi(x) = P(A = 1 | X = x)$, the conditional outcome density $p_{Y|X}(y | x)$, and the marginal covariate distribution $Q$. Let $m(x) = E_P[Y | X = x]$ denote the ‘outcome regression function’, such that the target estimand can be written as $\chi(P) = Q[m(X)]$. Thus, a plug-in estimator of $\chi$ takes the form $\chi(\hat{P}) = \hat{Q}[\hat{m}(X)]$ for estimates $\hat{Q}$ and $\hat{m}$ of $Q$ and $m$, respectively.

The efficient influence function of $\chi$ is

$$\hat{\chi}_P(Z) = \frac{A}{\pi(X)}\{Y - m(X)\} + m(X) - \chi(P), \quad (2)$$

which was first derived in Robins and Rotnitzky (1992). As a result, the one-step estimator can be easily constructed as

$$\hat{\chi}_{1-step} = \mathbb{P}_n \left[ \frac{A}{\hat{\pi}(X)} \{Y - \hat{m}(X)\} + \hat{m}(X) \right],$$

where we have avoided the need to estimate $Q$ but now require an estimate $\hat{\pi}$ of the propensity score. In Section 4, we will study this example in the Bayesian setting.

## 2.2 The one-step posterior correction

As discussed in the introduction, nonparametric smoothing in our prior may introduce a nonnegligible bias in the posterior on a low-dimensional functional that precludes a BvM theorem.
In some settings, it has been shown that this bias can be avoided by undersmoothing. But by doing so, the estimation of other aspects of the parameters can be negatively impacted and become sub-optimal (Castillo, 2012; Castillo & Rousseau, 2015; Ray & van der Vaart, 2020; Rousseau, 2016). Moreover, undersmoothing may require nontrivial modifications to existing implementations of Bayesian computation methods, as in Ray and van der Vaart (2020).

Instead of tuning the prior, we propose an approach that adjusts the posterior directly. Our starting point is an ordinary Bayesian model that is—as usual—optimized to recover the whole data-generating distribution $P$ rather than any particular functional. This means that the prior can be specified according to conventional guidelines, and the posterior $\Pi(\cdot | Z_{1:n})$ can be computed with standard methods.

**Definition 2** The one-step posterior of the functional $\chi$ given $Z_{1:n}$ is the pushforward measure of

$$ \tilde{\chi}: \mathcal{P} \times M(\mathcal{Z}) \to \mathbb{R} \\ \tilde{\chi}(P, \tilde{P}) = \chi(P) + \tilde{P}[\dot{\chi}_P]. $$

on the product measure

$$ (P, \tilde{P}) | Z_{1:n} \sim \Pi(\cdot | Z_{1:n}) \times \Pi_{BB}(\cdot | Z_{1:n}), $$

where $\Pi_{BB}$ denotes the Bayesian bootstrap posterior introduced in Section 2.1.

The structure of the corrected parameter $\tilde{\chi}$ takes inspiration from the one-step estimator described in the previous section. Informally speaking, the initial estimate $\tilde{P}$ has been replaced by the parameter $P$ drawn from the initial posterior, and the empirical measure $\mathbb{P}_n$ has been replaced by the Bayesian bootstrap parameter $\tilde{P}$. In our case, we are correcting an entire probability distribution rather than just a point estimator.

We can sample from the one-step posterior via the simple hierarchical scheme in Algorithm 1. First, the user draws samples $\{P^{(1)}, \dots, P^{(B)}\}$ of $P$ from their initial posterior $\Pi(\cdot | Z_{1:n})$ (e.g. via Markov chain Monte Carlo). Each posterior sample $P^{(b)}$ is then paired with an independent vector of uniform Dirichlet weights of length $n$; this arises from the representation

$$ \tilde{P} = \sum_{i=1}^{n} W_i \delta_{Z_i} $$

of the Bayesian bootstrap (e.g. Section 4.7 of S. Ghosal & van der Vaart, 2017), where $(W_1, \dots, W_n) \sim \text{Dir}(n; 1, \dots, 1)$. Finally, the draw of the estimand $\chi(P^{(b)})$ is summed with the weighted average of the efficient influence function at $P^{(b)}$.

Algorithm 1 is attractive for several reasons. First, the same set of posterior samples of $P$ can be retained to use for different functionals. Second, there is no need to modify any existing implementation for sampling from $\Pi(\cdot | Z_{1:n})$; the user can just take the output from the initial implementation and add the random correction terms. This extra step is extremely efficient as it only requires drawing independent sets of Dirichlet weights and averaging the efficient influence functions. Moreover, Steps 3–6 in Algorithm 1 can be easily parallelized to save further time.

To provide some geometric intuition for our method, we adapt the ideas from Fisher and Kennedy (2021). Consider the following mixture distributions constructed from an arbitrary $P \in \mathcal{P}$ (representing a draw from our initial posterior) and the true distribution $P_0$:

$$ P^{\epsilon} = (1 - \epsilon)P + \epsilon P_0 $$

for $0 \le \epsilon \le 1$. Within our model space $\mathcal{P}$, we could interpret $\{P^{\epsilon}\}_{\epsilon \in [0,1]}$ as a line segment linking $P$ and $P_0$ (these line segments are not guaranteed to lie within the model in general, but we will overlook this for the purpose of exposition). This induces a path $\{\chi(P^{\epsilon})\}_{\epsilon \in [0,1]}$ in the parameter space that reaches the true parameter value $\chi(P_0)$ at $\epsilon = 1$. Starting at $\epsilon = 0$, we might hope to get closer to the truth by constructing a good approximation to this path.
<caption>**Algorithm 1** One-step posterior sampling scheme</caption>

1 Input: number of posterior samples $B$;
2 Sample $\{P^{(1)}, \dots, P^{(B)}\}$ from the initial posterior $\Pi(\cdot | Z_{1:n})$;
3 For $b \leftarrow 1$ to $B$ do
4 &nbsp;&nbsp;&nbsp;&nbsp;Sample $(W_1^{(b)}, \dots, W_n^{(b)})$ from Dir($n$; 1, \dots, 1) independently of $P^{(b)}$;
5 &nbsp;&nbsp;&nbsp;&nbsp;Compute
6 end
7 Return $\{\tilde{\chi}^{(1)}, \dots, \tilde{\chi}^{(B)}\}$.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** One-step corrections as linear projections. The solid curve (—) is $\chi(P^\epsilon) - \chi(P_0)$, whose tangent at $\epsilon = 0$ is given by the dotted line (---). The dashed line (---) estimates the tangent with the empirical distribution, and the histogram (■) on the right corresponds to the distribution of the tangent drawn from the Bayesian bootstrap, representing the uncertainty in the estimand conditional on the starting point $P$.</figcaption></figure>

This idea is illustrated visually in Figure 1. The blue solid curve is equal to $\chi(P^\epsilon) - \chi(P_0)$, and its tangent at $\epsilon = 0$—given by the blue dotted line—has slope $P_0[\tilde{\chi}_P]$. The one-step estimator—taking $P$ as the ‘initial estimate’ of $P_0$—estimates the slope with $P_n[\tilde{\chi}_P]$, and the resulting projection onto $\epsilon = 1$ is shown by the green dashed line. With our posterior correction, the Bayesian bootstrap term $\tilde{P}[\tilde{\chi}_P]$ gives us a distribution on the slope rather than just a point estimate. As a result, we also get a distribution on the intercept at $\epsilon = 1$ after projecting along the tangent from $\epsilon = 0$, which is shown by the red histogram on the right. This intercept is the one-step corrected parameter; the histogram represents the uncertainty in $\tilde{\chi}$ conditional on the ‘posterior draw’ $P$. For the one-step posterior, this uncertainty is combined with the posterior uncertainty for $P$, which we recall to be conditionally independent of $\tilde{P}$ given the data.

The reader may wonder why the empirical measure is insufficient for our posterior correction. We can illustrate the reason by considering linear functionals, which form the simplest class of examples. Let $\chi(P) = P[g]$, where $g: \mathcal{Z} \to \mathbb{R}$ is a known bounded function. The efficient influence function is $\tilde{\chi}_P = g - P[g]$ (e.g. Castillo & Rousseau, 2015), such that

$$ \tilde{\chi} = P[g] + \tilde{P}[g - P[g]] = \tilde{P}[g]. $$
For this special case, the corrected parameter does not depend on the initial posterior at all: $\tilde{\chi}$ is simply the Bayesian bootstrap expectation of g, which we already know to satisfy the semiparametric BvM theorem (see, for example, Theorem 12.2 of S. Ghosal & van der Vaart, 2017 and the subsequent discussion). If $\tilde{P}$ is replaced by $P_n$ in the correction, then the parameter reduces to the one-step estimator $P_n[g]$, which has zero posterior uncertainty. Thus, the Bayesian bootstrap correction term serves not only to bias the initial posterior but also to obtain the (asymptotically) correct posterior shape. This statement is established in generality in Section 2.3.

In the preceding example of a linear functional, we highlight that it is in fact possible to use the Bayesian bootstrap for the initial posterior as linear functionals remain well-defined with respect to discrete distributions. This is also true for certain other functionals including quantiles and Z-functionals (Section 3.10.5.7 of van der Vaart & Wellner, 2023), which are defined as zeros of criterion functions. In these cases, the one-step posterior distribution of the target estimand is derived just from the Bayesian bootstrap, and the influence of the initial posterior $\Pi(\cdot | Z_{1:n})$ disappears. However, outside of these specific cases, the functionals of interest in semiparametric problems are usually defined with respect to restrictions on the data-generating distribution that require some degree of smoothness, in which case the Bayesian bootstrap is no longer feasible. For instance, we might require $\mathcal{P}$ to be dominated by the Lebesgue measure (e.g. Section 3), or we might impose positivity conditions as described in Example 1 (see also: Sections 4 and 5).

To summarize, the one-step posterior combines the uncertainty from two unknown quantities. First, we have the uncertainty in the functional $\chi(P)$, which arises from our initial posterior. Given a 'starting point' $P$, we also have uncertainty in $P_0[\dot{\chi}_P]$, which is the slope of the tangent to the curve connecting $\chi(P)$ to the truth $\chi(P_0)$. This latter source of uncertainty is modelled by replacing the unknown $P_0$ with the Bayesian bootstrap parameter $\tilde{P}$.

## 2.3 Theory

In this section, we provide conditions for the one-step posterior to satisfy the semiparametric BvM theorem; that is,

$$d_{BL}(\mathcal{L}_{\Pi \times \Pi_{BB}}(\sqrt{n}(\tilde{\chi} - \hat{\chi}_n) | Z_{1:n}), \mathcal{N}(0, \|\dot{\chi}_{P_0}\|_{P_0}^2)) \xrightarrow{P_0} 0, \quad (3)$$

where $\mathcal{L}_{\Pi \times \Pi_{BB}}(\sqrt{n}(\tilde{\chi} - \hat{\chi}_n) | Z_{1:n})$ denotes the posterior law of $\sqrt{n}(\tilde{\chi} - \hat{\chi}_n)$ and $\hat{\chi}_n$ is the asymptotically efficient sequence $\hat{\chi}_n = \chi(P_0) + P_n[\dot{\chi}_{P_0}]$. As discussed in Section 2.1, this justifies the use of central credible sets as confidence regions.

**Assumption 1** There exists a sequence of measurable subsets $(H_n)_{n \in \mathbb{N}}$ of $\mathcal{P}$ satisfying $\Pi(P \in H_n | Z_{1:n}) \xrightarrow{P_0} 1$ with

(a) (No second-order bias)

$$\sup_{P \in H_n} |\sqrt{n}r_2(P_0, P)| = \sup_{P \in H_n} |\sqrt{n}\{\chi(P_0) - \chi(P) - P_0[\dot{\chi}_P]\}| \to 0.$$

(b) ($L_2$-convergence) $\sup_{P \in H_n} \|\dot{\chi}_P - \hat{\chi}_{P_0}\|_{P_0} \to 0$

(c) (Donsker class) The sets $\{\dot{\chi}_P : P \in H_n\}$ are eventually contained in a fixed $P_0$-Donsker class (van der Vaart & Wellner, 1996).

We also provide an alternative to Assumption 1(c) below: $1(c^*)(i)$ weakens the Donsker condition at the cost of adding $1(c^*)(ii)$, which requires envelope function moment bounds. This alternative is also more appropriate if the sieve $\{H_n\}$ is not nested (i.e. we do not have ... $\supset H_{k-1} \supset H_k \supset H_{k+1} \supset ...$), which would likely preclude the possibility of finding a fixed Donsker class that satisfies Assumption 1(c).
**Assumption 1** (c*) (i) (Convergence of $\dot{\chi}_p$ under the empirical process)

$$ \sup_{P \in H_n} | \mathbb{G}_n[\dot{\chi}_P - \dot{\chi}_{P_0}] | \xrightarrow{P_0} 0. $$

This statement is interpreted in terms of outer probability (van der Vaart & Wellner, 1996) if the supremum is not measurable.

(ii) (Bounding of envelope functions) The sets $\{\dot{\chi}_p: P \in H_n\}$ have envelope functions $G_n$ (i.e. $|\dot{\chi}_p(z)| \le G_n(z)$ for all $P \in H_n$ and all $z \in \mathcal{Z}$) satisfying

$$ \lim_{C \to \infty} \limsup_{n \to \infty} P_0 G_n^2 1_{G_n^2 > C} = 0, \quad P_0 G_n^4 = o(n). $$

**Theorem 1** Under Assumption 1, the one-step posterior satisfies the semiparametric BvM theorem so that (3) holds true. The result (3) also holds if Assumption 1(c) is replaced by Assumption 1(c*).

The proof of Theorem 1 can be found in online supplementary material, Section S8.1. The structure of Assumption 1 mirrors the standard assumptions for one-step estimators to attain asymptotic efficiency (see online supplementary material, Section S1.2). To satisfy Assumptions 1(a) and (b), the sets $H_n$ are typically included in shrinking neighbourhoods of $P_0$. In light of the earlier discussion about the quadratic nature of $r_2$, this term can be neglected when these neighbourhoods shrink quickly enough. The empirical process conditions in Assumptions 1(c) and (c*) restrict the complexity of the posterior so that we can obtain uniform convergence results.

To provide some context for Assumption 1, we compare with the conditions for Theorem 2.1 in Castillo and Rousseau (2015), which is a semiparametric BvM result for functionals. As one might expect, the second-order remainder must be negligible in both cases [Assumption 1(a)] so that the functional is well-approximated by its linear expansion. Assumptions 1(b), (c), and (c*) are not explicitly required by Castillo and Rousseau (2015), but similar conditions are generally used to establish the asymptotically normal expansion of the likelihood [see (4.10) of Castillo & Rousseau, 2015, Theorem 1 of Ray & van der Vaart, 2020 and p. 373–374 of S. Ghosal & van der Vaart, 2017]. Analogous conditions can also be found in the classical setting for semiparametric likelihood methods (Cheng & Kosorok, 2008b; Murphy & van der Vaart, 2000).

Crucially, our result avoids the restrictive ‘change-of-measure’ condition in (2.14) of Castillo and Rousseau (2015) [see also (3.6) in Ray & van der Vaart, 2020]. This depends delicately on the interplay between the model and the prior, which can fail if the prior is oversmooth or adaptive. From a technical viewpoint, the condition is also challenging to verify for complex priors.

The crux of our proof lies in the decomposition

$$ \sqrt{n}(\tilde{\chi} - \hat{\chi}_n) = \underbrace{\sqrt{n}(\tilde{P} - \mathbb{P}_n)[\dot{\chi}_P]}_{(I)} + \underbrace{\mathbb{G}_n[\dot{\chi}_P - \dot{\chi}_{P_0}]}_{R_n(P)} - \sqrt{n}r_2(P_0, P). \quad (4) $$

Roughly speaking, the contribution of the remainder term $R_n(P)$ is negligible to first-order under the assumptions of the theorem, and the asymptotic behaviour of the posterior law is driven by the leading term (I). To study term (I), the conditional independence of $P$ and $\tilde{P}$ given $Z_{1:n}$ is crucial because we wish to utilize the desirable asymptotic properties of the Bayesian bootstrap, and these properties must continue to hold even if we condition on $P$. This conditional independence resembles the use of sample-splitting in classical semiparametric statistics (Chernozhukov et al., 2018; Hines et al., 2022) but serves a different purpose. Sample-splitting and cross-fitting are employed to weaken the requisite empirical process conditions for efficient inference (see online supplementary material, Section S1.2 for further details). This is closely linked to the example of Ray and van der Vaart (2020) that we study in Section 4.3, where the propensity score is estimated on an auxiliary dataset. On the other hand, $P$ and $\tilde{P}$ are independent conditional on the same dataset, so we do not obtain a similar relaxation of empirical process assumptions. A possible analogous approach would be to fit the initial and Bayesian bootstrap posteriors on separate splits
of the data, which would raise the possibility of relaxing Assumptions 1(c) and (c*). We leave the exploration of this idea to future work.

Interestingly—as noted in the proof of Theorem 1—our approach can be used on any distribution $\hat{\Pi}$ that satisfies the assumptions, so $\hat{\Pi}$ does not need to be a proper posterior distribution. More precisely, let $\hat{\Pi}$ be any distribution on $\mathcal{P}$ such that $\hat{\Pi}(H_n) = 1 + o_{P_0}(1)$ for $H_n$ satisfying Assumption 1. Then

$$d_{BL}(\mathcal{L}_{\hat{\Pi} \times \Pi_{BB}}(\sqrt{n}(\tilde{\chi} - \hat{\chi}_n) | Z_{1:n}), \mathcal{N}(0, \|\dot{\chi}_{P_0}\|_{P_0}^2)) \xrightarrow{P_0} 0. \quad (5)$$

The result is particularly interesting for variational Bayes posteriors (Jordan et al., 1999) because apart from specific setups (e.g. when the posterior satisfies a BvM and the variational family is Gaussian), the variational Bayes measures of uncertainty have at best no justification and at worst are misleading (e.g. mean field variational Bayes). There are now several results on variational posterior contraction rates (Alquier & Ridgway, 2019; Zhang & Gao, 2020) to help verify the assumptions. For example, using Lemma 15 of Sulem et al. (2023), if $\Pi(H_n^c | Z_{1:n})$ goes to zero sufficiently quickly on sets of $P_0^n$-probability going to 1, then the variational Bayes posterior satisfies $\hat{\Pi}(H_n) = 1 + o_{P_0}(1)$ under the same assumptions as in Zhang and Gao (2020).

Another interesting example arises if we set a Dirac mass posterior $\hat{\Pi}(\cdot) = \delta_p(\cdot)$ for a chosen point estimator $\hat{P}$, e.g. the posterior mean of $P$ for an initial posterior. The one-step posterior mean is then equal to the one-step estimator $\chi(\hat{P}) + P_n[\dot{\chi}_P]$, and we derive the one-step posterior by ‘Bayesian bootstrapping’ the one-step estimator with $\hat{P}$ fixed, akin to the percentile bootstrap (Efron & Hastie, 2016). We compare this with the classical approach for constructing Wald confidence intervals for the one-step estimator (e.g. Hines et al., 2022), which also keeps $\hat{P}$ fixed. Here, the variance of the one-step estimator is estimated with the sample variance of the estimated efficient influence function, $n^{-1}P_n[(\dot{\chi}_P - P_n[\dot{\chi}_P])^2]$.

The sample variance method outlined above is based on a decomposition of the one-step estimator similar to (4) (see online supplementary material, Section S1.2 or Hines et al., 2022). Its validity depends in particular on remainder terms being higher-order so that they can be ignored in the variance calculations; notably, the uncertainty associated with $\hat{P}$ is ignored. In the small sample regime, however, the remainder terms can be nonnegligible, in which case the interval estimates are at risk of being anticonservative. This is a well-known phenomenon in semiparametric statistics (e.g. Gruber & van der Laan, 2019; van der Laan et al., 2021). Owing to the complex behaviour of flexible black-box algorithms, it is generally impossible to remedy this by bootstrapping the point estimator $\hat{P}$ (see p. xviii and Chapter 28 of van der Laan & Rose, 2018). We argue that our one-step posterior framework is advantageous for small sample inference because the uncertainty on $P$—both in $\dot{\chi}_P$ and in $R_n(P)$ in (4)—is accounted for by the posterior distribution on $P$. With our simulation results in Section 6.2, we verify that the two plug-in approaches—the Dirac mass one-step posterior and frequentist sample variance estimation—do indeed match each other very closely. More importantly, the results illustrate our point that the plug-in approaches are susceptible to considerable undercoverage in finite sample inference, while the one-step posterior can incorporate the additional uncertainty in the remainder terms and perform significantly better.

A formal theoretical justification of our arguments would potentially require a higher-order analysis. There are some precedents for this type of theory for semiparametric Bayes-like methods (Cheng & Kosorok, 2008a, 2008b); we leave this as a future research direction. In any case, we draw a comparison with the role of the prior distribution in parametric Bayes. The parametric BvM theorem asserts that the posterior derived from a parametric model converges to the sampling distribution of the maximum likelihood estimator, and the prior only has a higher-order contribution (e.g. Chapter 10 of van der Vaart, 1998). However, it is well-known that the prior can nonetheless play a key role in finite-sample inference (p. 125, Point 4 of J. Berger, 1985). Our philosophy aligns with the ‘Calibrated Bayes’ perspective of Rubin (1984) and Little (2011): a calibrated Bayesian procedure offers the potential for improved finite-sample inference but allows the data to dominate as the sample size increases.

It is possible to replace the Bayesian bootstrap with proper Dirichlet process posteriors but this would necessitate additional assumptions on the prior base measure, so we will not pursue this further. We also expect that it is viable to replace the Bayesian bootstrap with Efron’s
nonparametric bootstrap, which has similar first-order asymptotic properties (e.g. Chapter 3.7 of van der Vaart & Wellner, 2023). The nonparametric bootstrap is used in a related approach by Antonelli et al. (2022), which shares our motivation of improving small sample inference using posterior uncertainty. The proposal of Antonelli et al. (2022) is based on the average treatment effect estimation problem, which is a direct generalization of our mean outcome example in Example 1 and Section 4. Their primary objective is to estimate the frequentist variance of a point estimator that coincides with our one-step posterior mean, and they rely on asymptotic normality of this estimator to construct confidence intervals. Their variance estimator is structurally different from the variance of the one-step posterior distribution, as explained in online supplementary material, Section S4, where further discussion on the differences between both approaches is provided. Note in particular that their approach targets a single estimand. We also provide numerical comparisons between both approaches in Section 6.2, where the confidence regions based on Antonelli et al. (2022) have a lower frequentist coverage, particularly for smaller sample sizes.

As we will see with our examples, the efficient influence function often admits a natural re-presentation of the form

$$ \dot{\chi}_P = \psi_P - P[\psi_P], $$

where $\psi_P: \mathcal{Z} \to \mathbb{R}$ is a measurable function. Consider Example 1 for concreteness, where we can define

$$ \psi_P(Z) = \frac{A}{\pi(X)} \{Y - m(X)\} + m(X), \quad (6) $$

which satisfies $P[\psi_P] = \chi(P)$. Theorem 1 still holds if $\dot{\chi}_P$ and $\dot{\chi}_{P_0}$ are replaced by $\psi_P$ and $\psi_{P_0}$, respectively, in Assumptions 1(b), (c), and (c*). This has the potential to simplify the verification of the assumptions.

### 3 The integrated squared density

Suppose that we observe independent and identically distributed data $Z_{1:n} = (Z_1, \dots, Z_n)$ from a distribution $P_0$ on $\mathbb{R}$ that is absolutely continuous with respect to Lebesgue measure with a square-integrable density. We make slight changes in the notation by identifying distributions with their Lebesgue densities, e.g. $P_0$ is now identified by the true density $f_0$. Consider the problem of estimating $\chi(f) = \int f(z)^2 dz$. There is a large literature on the estimation of $\chi(f)$ (e.g. Bickel & Ritov, 1988; Birgé & Massart, 1995; Giné & Nickl, 2008; Laurent, 1996; Newey et al., 1998) and more generally, on power integrals $\int f(z)^q dz$ for integers $q \ge 2$. The integrated squared density plays an important role in nonparametric rank statistics (p. 266–267 of Prakasa Rao, 1983).

In Castillo and Rousseau (2015), the authors obtained a BvM theorem for this functional with random histogram priors, which are much simpler than common priors used in density estimation such as infinite mixture models. Under regularity assumptions on the true density $f_0$, we prove in this section that the one-step posterior distribution satisfies the semiparametric BvM theorem for Dirichlet process Gaussian location mixtures. The prior on $f$ is defined by

$$ f_{G,\sigma}(z) = \int_{\mathbb{R}} \varphi_{\sigma}(z - \mu) dG(\mu), \quad G \sim DP(MG_0), \quad M > 0, \quad (7) $$

where $\varphi_{\sigma}(z - \mu)$ is the Gaussian density centred at $\mu$ with standard deviation $\sigma$, and $G_0$ is a distribution on $\mathbb{R}$ with density $g_0(x)$ that is positive and continuous and satisfies $g_0(x) \lesssim e^{-c_1|x|^\tau}$ for some $\tau > 0$. For instance, we can choose $G_0$ to be a Gaussian distribution. As in Kruijer et al. (2010), Shen et al. (2013), and Scricciolo (2014), we consider an inverse Gamma prior on $\sigma: \sigma \sim IG(a, b)$.

The efficient influence function for $\chi$ is $\dot{\chi}_f(Z) = 2\{f(Z) - \chi(f)\}$ (Bickel & Ritov, 1988), so the one-step corrected parameter is

$$ \tilde{\chi} = 2\tilde{P}[f(Z)] - \chi(f). $$
It is straightforward to derive the second-order bias:

$$r_2(f_0, f) = \|f - f_0\|_2^2,$$

where $\| \cdot \|_2$ denotes the $L_2$-norm with respect to Lebesgue measure. We highlight the quadratic nature of this term; accordingly, we require the $L_2$ contraction rate for $f$ to be faster than $n^{-1/4}$.

We consider the following smoothness assumption on $f_0$ introduced by Shen et al. (2013).

**Definition 3** A density $f$ belongs to the class of locally $\alpha$-Hölder functions $C(\alpha, L, \tau_0)$ with function $L$ and $\tau_0 \ge 0$ if $p$ has $k_\alpha = [\alpha] - 1$ derivatives with

$$|(D^{k_\alpha} f)(z+h) - (D^{k_\alpha} f)(z)| \le L(z)e^{\tau_0 h^2} |h|^{\alpha-k_\alpha}, \quad \forall z, h \in \mathbb{R}.$$

**Theorem 2** Suppose that the true density $f_0$ belongs to a class $C(\alpha, L, \tau_0)$ for $\alpha > 1$. Assume further that there exists $\epsilon > 0$ such that for all $k \le k_\alpha$

$$P_0\left[\left(\frac{|D^k f_0|}{f_0}\right)^{(2\alpha+\epsilon)/k}\right] < \infty, \quad P_0\left[\left(\frac{L}{f_0}\right)^{2+\epsilon/\alpha}\right] < \infty$$

and $f_0$ is monotone nonincreasing (resp. nondecreasing) for $z$ large enough (resp. $-z$ large enough) and there exist $c_1, c_2 > 0$ such that

$$f_0(z) \le c_1 e^{-c_2|z|^{\tau_0}}.$$

If we specify a Dirichlet process Gaussian location mixture prior for $f$ as defined above, then the one-step posterior distribution of the integrated squared density satisfies the semiparametric BvM theorem.

The proof of Theorem 2 can be found online supplementary material, Section S8.2. From Shen et al. (2013), under the above conditions the posterior distribution concentrates in Hellinger or $L_1$ distance at the rate $n^{-\alpha/(2\alpha+1)}(\log n)^r$ for some $r > 0$, which is minimax adaptive up to a $\log n$ term for all $\alpha > 0$ (Maugis-Rabusseau & Michel, 2013). The condition $\alpha > 1$ is used to verify part (c*) of Assumption 1. It can be weakened to $\alpha > 1/2$ by considering the truncated prior: $f \sim \Pi_T(df) \propto \Pi(df) 1_{\|f\|_\infty < K}$ for some arbitrarily large but fixed constant $K$.

Note that there are different definitions of the class $C(\alpha, L, \tau_0)$ in the literature, and all these variants can be used in Theorem 2. As in Donnet et al. (2018), this condition can be replaced by assuming the existence of a finite mixture of Gaussian densities $f_{G^*,\sigma}$ with variance $\sigma^2 = o(1)$ with at most $N = O(\sigma^{-1} \log 1/\sigma)$ support points such that the Kullback-Leibler divergence between $f_0$ and $f_{G^*,\sigma}$ is of order $O(\sigma^{2\alpha})$.

# 4 The mean of an outcome that is missing-at-random

## 4.1 General methodology

In this section, we revisit the missing data example in Example 1. Recall that the data take the form $Z = (X, A, AY)$, where $X$ is a vector of covariates, $A$ is a binary missingness indicator, $Y$ is the outcome variable of interest that is unobserved when $A = 0$, and the target estimand is

$$\chi(P) = E_P[E_P[Y | A = 1, X]].$$

The density for a single observation $z = (x, a, ay)$ can be written as

$$p(z) = \pi(x)^a (1 - \pi(x))^{1-a} p_{Y|X}(y | x)^a q(x),$$

which factorizes over the parameters $(\pi, p_{Y|X}, q)$. Consequently, if the three parameters are given independent priors, they will remain independent given $Z_{1:n}$ in the posterior. An important
implication is that the model for $\pi$ will not feature in the marginal posterior of $\chi$ since the functional is only defined by the other two parameters. This well-known phenomenon is discussed in Robins and Ritov (1997), Robins et al. (2015), and Ray and van der Vaart (2020).

The efficient influence function of $\chi$ is given by (2) so that the one-step corrected parameter

$$\tilde{\chi}(P, \tilde{P}) = \tilde{P} \left[ \frac{A}{\pi(X)} \{Y - m(X)\} + m(X) \right]$$

depends on $P$ only through $\pi$ and $m$. Since $Q$ is independent of $\pi$ and $m$ in the initial posterior, the model specification for the marginal covariate distribution does not affect the corrected parameter.

In the following, we state assumptions for the one-step posterior to satisfy the semiparametric BvM theorem, specializing the conditions in Assumption 1. Let $\pi_0$ and $m_0$, respectively, denote the true values of $\pi$ and $m$ associated with $P_0$.

**Assumption 2** There exists a sequence of measurable subsets $(H_n)_{n \in \mathbb{N}}$ of $P$ satisfying $\Pi(P \in H_n | Z_{1:n}) \xrightarrow{n} 1$ such that

(a) (L2-convergence of $\pi$ and $m$) There exist numbers $\rho_n, \epsilon_n \to 0$ such that

$$\sup_{P \in H_n} \left\| \frac{1}{\pi} - \frac{1}{\pi_0} \right\|_{P_0} \le \rho_n, \quad \text{and} \quad \sup_{P \in H_n} \|m - m_0\|_{P_0} \le \epsilon_n$$

$$\text{and } \sqrt{n}\rho_n\epsilon_n \to 0.$$

(b) (Uniform bounding) For all sufficiently large $n$, there exist fixed constants $\delta, C > 0$ such that for all $P \in H_n$,

$$\left\| \frac{1}{\pi} \right\|_{\infty} + \left\| \frac{1}{\pi_0} \right\|_{\infty} < \frac{1}{\delta}, \quad \text{and} \quad P_0(|Y - m(X)| < C) = 1.$$

(c) (Donsker class) The sequences of sets $\{\pi: P \in H_n\}$ and $\{m: P \in H_n\}$ are both eventually contained in fixed $P_0$-Donsker classes.

**Theorem 3** Under Assumption 2, the one-step posterior for the outcome mean satisfies the semiparametric BvM theorem.

The proof of Theorem 3 can be found in online supplementary material, Section S8.3. In online supplementary material, Section S5.1, we provide an example of verifying Assumption 2 by studying the binary nonparametric regression prior from Ghosal and van der Vaart (2007). Assumption 2(a) states that $1/\pi$ and $m$ must both converge uniformly on $H_n$ to their respective truths in $L_2$, and their combined rate of convergence must be faster than $1/\sqrt{n}$. To connect this to Assumption 1(a), we note that the second-order bias takes the form

$$r_2(P_0, P) = P_0 \left[ \pi_0(X) \left( \frac{1}{\pi_0(X)} - \frac{1}{\pi(X)} \right) \{m_0(X) - m(X)\} \right] \le \left\| \frac{1}{\pi_0} - \frac{1}{\pi} \right\|_{P_0} \|m_0 - m\|_{P_0}, \quad (8)$$

where the upper bound is due to Cauchy-Schwarz. Thus, Assumption 2(a) immediately implies Assumption 1(a). The cross-term structure of the second-order bias enables the rate of convergence of one parameter to compensate for the other; this is a property known as ‘rate double robustness’ (Rotnitzky et al., 2021). However, attaining rate double robustness is complicated by the requisite empirical process condition (c) in Assumption 2; we discuss this further in Section 4.2.

The uniform bounds in Assumption 2(b) help to ensure that the $L_2$-convergence of $\chi_P$ [i.e. Assumption 1(b)] follow from Assumption 2(a). The assumption that the propensity score is bounded away from zero is commonly referred to as ‘positivity’ and is routine for any estimation approach that involves inverse probability weighting.
## 4.2 Double robustness

In Section 4.1, we briefly discussed the ‘rate double robustness’ structure of Assumption 2(a), highlighting that the product of the $L_2(P_0)$-contraction rates for $\pi$ and $m$ must be faster than $n^{-1/2}$. For example, the condition would be satisfied provided that the posteriors for $\pi$ and $m$ both contract faster than $n^{-1/4}$. The more interesting setting arises when one of $(\pi, m)$ contracts slower than $n^{-1/4}$ but the other parameter compensates by contracting faster. In this section, we provide an explicit example of attaining the BvM theorem in this setting.

Our mean outcome example was studied by Ray and van der Vaart (2020) under the constraint that the outcome variables are binary and $X \in [0, 1]^d$ ($d \ge 1$). We will work in their set-up with $d = 1$ and compare with their results later on. The outcome regression function is modelled by $m(x) = \Psi(w(x))$, where $\Psi(t) = 1/(1 + e^{-t})$ is the expit function, and the parameter $w$ follows a Riemann-Liouville process prior released at zero of regularity $\bar{\beta} > 0$, defined by

$$w(x) = \sum_{k=0}^{|\bar{\beta}|+1} g_k x^k + \int_0^x (x-s)^{\bar{\beta}-1/2} dB_s, \quad x \in [0, 1]. \quad (9)$$

In the above display, $|\bar{\beta}| + 1$ is the smallest integer that is strictly larger than $\bar{\beta}$, the {$g_k$} are i.i.d. $N(0, 1)$ variables, and $B$ is an independent Brownian motion (see Chapter 11 of S. Ghosal & van der Vaart, 2017).

For the propensity score model, we specify a simpler sequence of priors. Let ($I_j, j \le J_a$) be a regular partition of $[0, 1]$, so that $|I_j| = J_a^{-1}, J_a \ge 2$. Assume that $\pi_j \stackrel{\text{iid}}{\sim} g$, where $g$ is a positive and continuous density on $[0, 1]$, and define

$$\pi(x) = \sum_{j=1}^{J_a} 1_{I_j}(x)\pi_j, \quad \pi_{0,j}(x) = \sum_{j=1}^{J_a} 1_{I_j}(x)\pi_{0j},$$

where $\pi_{0j} = p_{0j}/Q_0(I_j)$ with $p_{0j} = \int_{I_j} \pi_0(x)Q_0(dx)$ for all $j \le J_a$. As described in Section 4.1, the priors on $(m, \pi, Q)$ are enforced to be mutually independent, so the model for $Q$ drops out of the analysis for the one-step posterior. We then have the following proposition.

**Proposition 1** Consider the prior defined above and assume that $\pi_0 \in \mathcal{C}^a([0, 1])$ with $0 < \alpha \le 1$ and $m_0 \in \mathcal{C}^\beta([0, 1])$, where $\mathcal{C}^s([0, 1])$ is the Hölder space of order $s$ on $[0, 1]$ (see Chapter 4 of Giné & Nickl, 2016). Assume also that there exists $c_0 > 0$ such that $c_0 \le \pi_0 \le 1 - c_0$. Then suppose $J_a = o(\sqrt{n/\log n})$, and $\rho_n$ and $\epsilon_n$ in Assumption 2 are given by

$$\rho_n = \left( J_a^{-\alpha} \vee \sqrt{\frac{J_a \log n}{n}} \right) \quad \text{and} \quad \epsilon_n = n^{-\frac{\beta_n \bar{\beta}}{2\beta+1}} \log n. \quad (10)$$

Moreover, if $\beta \wedge \bar{\beta} > 1/2$, $\sqrt{n}\rho_n\epsilon_n \to 0$ and if for some $t > 0$

$$J_a^3 = O\left( n^{\frac{2\bar{\beta}}{2\beta+1} + \frac{2(\bar{\beta} - \bar{\beta}_n)\beta(2\bar{\beta} - 1)}{2\beta + 1} - t \right), \quad (11)$$

then the one-step posterior for the outcome mean satisfies the semiparametric BvM theorem.

The proof of Proposition 1 is provided in online supplementary material, Section S8.4. It follows a slightly different strategy to Theorem 3 because it replaces the Donsker class condition in Assumption 2(c) with an adapted version of Assumption 1(c*). This involves controlling $|\mathbb{G}_n[\psi_p - \psi_{P_0}]|$ uniformly over an appropriate sequence of subsets, where $\psi_p$ was defined in (6).

We compare our result with Corollary 1 from Ray and van der Vaart (2021), which studied a fully Bayesian analysis under the same set-up with the same Riemann-Liouville prior on $w$ and an independent Dirichlet process prior on $Q$ (the propensity score drops out of the analysis because it
is given an independent prior, as explained in Section 4.1). They required $\alpha, \beta > 1/2$ and $1/2 < \bar{\beta} < \alpha + \beta - 1/2$ to attain the semiparametric BvM theorem. In contrast, we see from (10) that the contraction rate $\rho_n$ of $\pi$ is slower than $J_a^{-\alpha}$. In particular, if we set $J_a = o(n^\delta)$ for $\delta < 1/2$, then $n^{\delta a}\rho_n \to \infty$. Thus, we can have $\alpha < 1/2$ and $n^{1/4}\rho_n \to \infty$ as long as $\beta \wedge \bar{\beta}$ is sufficiently large to compensate. This yields one half of rate double robustness.

To better understand Proposition 1, the key reason behind this one sided double robustness is that when $J_a \lesssim n^{-\delta}$ with $\delta$ small (i.e. $J_a$ increases slowly to infinity), the complexity of the prior is not so large but it leads to a larger bias, in the sense that there exists $\pi_{0,J_a}(x) = \sum_{j=1}^{J_a} 1_{I_j}(x)\pi_{0j}$ such that the posterior concentrates at the fast rate $\sqrt{J_a \log n/n}$ around $\pi_{0,J_a}$. Then Assumption 1(c*) is verified by controlling $\mathbb{G}_n[1/\pi - 1/\pi_{0J_a}]$ using the fact that $J_a$ is not too large. When the smoothness $\alpha$ of $\pi_0$ is small, the posterior contraction rate around $\pi_0$ is dominated by the bias term of order $J_a^{-\alpha}$, which does not affect the control of Assumption 1(c*). We conjecture that a similar phenomenon takes place with the Riemann-Liouville process if $\beta > 1/2$ while $\beta < \bar{\beta}$, but this would require a much more technical proof and is left for future investigation. This would imply in particular a full double robustness behaviour under Gaussian process priors.

## 4.3 Comparison with prior augmentations in the least favourable direction

In this section, we compare our correction method with a proposal by Ray and van der Vaart (2020) that has similar motivations. As described in Section 4.2, Ray and van der Vaart (2020) studied the same outcome mean example under the constraint that the outcome variables are binary. They proposed constructing a preliminary estimator $\hat{\pi}_n$ of the propensity score based on a separate dataset that is independent of $Z_{1:n}$ and modelling the outcome regression function by

$$m(x) = \Psi\left\{w(x) + \frac{\lambda}{\hat{\pi}_n(x)}\right\}, \qquad (12)$$

where $\Psi(t) = 1/(1 + e^{-t})$ is the expit function, and $(w, \lambda)$ are parameters.

The prior on $w(x)$ should be specified such that the resulting prior on $\Psi\{w(x)\}$ captures the user's prior beliefs on $m(x)$. In contrast, the augmentation term $\lambda/\hat{\pi}_n(x)$ plays a technical role, with the goal of satisfying the change-of-measure condition for the semiparametric BvM theorem as discussed in Section 2.2. The fluctuation parameter $\lambda$ follows a $N(0, \sigma_n^2)$ prior independently of $w(x)$, where the variance $\sigma_n^2$ is allowed to vary with $n$. For the marginal covariate distribution, Ray and van der Vaart (2020) specified a Dirichlet process prior independent of $(w(x), \lambda)$. We will assume throughout that $\pi_0$ and $\hat{\pi}_n$ satisfy the following assumption.

**Assumption 3** The true propensity score satisfies $\pi_0 > \delta$ for some $\delta > 0$, and the preliminary estimators $\hat{\pi}_n$ satisfy $\|1/\hat{\pi}_n\|_{\infty} = O_{P_0}(1)$ and there exists $\rho_n \to 0$ such that

$$\left\| \frac{1}{\hat{\pi}_n} - \frac{1}{\pi_0} \right\|_{P_0} = O_{P_0}(\rho_n),$$

**Theorem 4** (Theorem 2 of Ray & van der Vaart, 2020). Assume that there exist measurable sets $\mathcal{H}_n$ of functions satisfying, for every $t \in \mathbb{R}$ and some numbers $u_n, \varepsilon_n \to 0$,

$$\Pi(\lambda: |\lambda| \le u_n \sigma_n^2 \sqrt{n} | Z_{1:n}) \xrightarrow{P_0} 1, \qquad (13)$$

$$\Pi((w, \lambda): w + (\lambda + tn^{-1/2})\hat{\pi}_n^{-1} \in \mathcal{H}_n | Z_{1:n}) \xrightarrow{P_0} 1, \qquad (14)$$

$$\sup_{m=\Psi(\eta): \eta \in \mathcal{H}_n} \|m - m_0\|_{P_0} \le \varepsilon_n, \qquad (15)$$

$$\sup_{m=\Psi(\eta): \eta \in \mathcal{H}_n} |\mathbb{G}_n[m - m_0]| \xrightarrow{P_0} 0. \qquad (16)$$

Under Assumption 3, if $n\sigma_n^2 \to \infty$ and $\sqrt{n}\rho_n\varepsilon_n \to 0$, then the marginal posterior of $\chi$ satisfies the semiparametric BvM theorem.
The critical aspect of the above result is that it avoids the change-of-measure condition in equation (3.6) of Theorem 1 of Ray and van der Vaart (2020) for establishing the semiparametric BvM theorem without the prior augmentation. We now compare this approach with the one-step posterior in this set-up. As discussed by Ray and van der Vaart (2020), the preliminary estimator $\hat{\pi}_n$ could be interpreted as a degenerate prior on $\pi$. This interpretation immediately suggests a corresponding one-step posterior by defining

$$ \tilde{\chi} = \tilde{P} \left[ \frac{A}{\hat{\pi}_n(X)} \{Y - m(X)\} + m(X) \right]. \quad (17) $$

To facilitate comparisons with the method of Ray and van der Vaart (2020), we will continue to describe the assumptions in terms of the parameter $w(x)$, where we have $m = \Psi(w)$ in this case.

**Theorem 5** Assume that there exist measurable sets $\tilde{H}_n$ of functions satisfying, for some numbers $\varepsilon_n \to 0$,

$$ \Pi(w \in \tilde{H}_n | Z_{1:n}) \xrightarrow{P_0} 1, \quad (18) $$

$$ \sup_{m=\Psi(w): w \in \tilde{H}_n} \|m - m_0\|_{P_0} \le \varepsilon_n, \quad (19) $$

$$ \sup_{m=\Psi(w): w \in \tilde{H}_n} |G_n[m - m_0]| \xrightarrow{P_0} 0. \quad (20) $$

Under Assumption 3, if $\sqrt{n}\rho_n\varepsilon_n \to 0$, then the one-step posterior satisfies the semiparametric BvM theorem.

The proof of Theorem 5 can be found in online supplementary material, Section S8.5. A direct comparison of the assumptions for Theorems 4 and 5 is complicated slightly by the different specifications for the outcome regression function, which means that the sets $H_n$ and $\tilde{H}_n$ may not match. Nevertheless, we should keep in mind that the augmentation in (12) is motivated purely by improving the inference for $\chi$ rather than the recovery of $m$. It is likely, therefore, that a prior on $w$ that satisfies the assumptions (13)–(16) (in combination with some sequence of priors for $\lambda$) will also satisfy the assumptions (18)–(20) for the unaugmented regression model. (In any case, it will be simpler to check (18)–(20) due to the absence of the augmentation term.) We verify this claim below for a Gaussian process prior studied by Ray and van der Vaart (2020).

Following the set-up in Ray and van der Vaart (2020), we take the covariate space to be $[0, 1]$. The parameter $w$ is given the Riemann–Liouville process prior defined in (9). For covariate spaces equal to $[0, 1]^d$ for integer $d \ge 1$, the subsequent corollary can be easily extended to the Gaussian series priors considered by Ray and van der Vaart (2020).

**Corollary 1** Suppose that $\pi_0^{-1} \in \mathbb{C}^a([0, 1])$, $m_0 \in \mathbb{C}^\beta([0, 1])$ for $a, \beta > 0$. Suppose further that we specify our one-step posterior using (17) and the Riemann–Liouville prior for $m$ defined by (9). Under Assumption 3, if $\beta \wedge \tilde{\beta} > 1/2$ and $\sqrt{n}\rho_n\varepsilon_n \to 0$ for

$$ \varepsilon_n = n^{-\frac{\beta n \tilde{\beta}}{2\beta+1}} \log n, $$

then the one-step posterior satisfies the semiparametric BvM theorem.

The corresponding result for the prior augmentation method in Corollary 2 of Ray and van der Vaart (2020) requires the following additional condition:

$$ \sigma_n = O(1) \quad \text{and} \quad \varepsilon_n \ll \sigma_n. $$

For finite-sample inference, the need to tune $\sigma_n$ is a substantial practical drawback of the prior augmentation approach relative to the one-step posterior. This is because the prior for $\lambda$ will likely
have a significant impact on the posterior for the target estimand, yet the augmentation term only plays a technical role, so there will be no subjective or substantive guidance to select $\sigma_n$.

From a computational perspective, the prior augmentation is also more difficult to implement than the one-step posterior correction—it is generally the case that existing software for posterior sampling will have to be rewritten in order to accommodate the new linear term. Furthermore, as Ray and van der Vaart (2020) wrote in their discussion, the augmented prior ‘will not perform better for any other functional’. Therefore, if multiple functionals are of interest, the user would have the time-consuming task of implementing distinct models for each functional. As we have mentioned previously, our approach handles this situation very efficiently because we can retain a single set of posterior samples from the initial model, with which we can compute all requisite corrected posteriors.

We recall that the theory provided by Ray and van der Vaart (2020) for this prior augmentation relies on an independent, auxiliary dataset to construct the propensity score estimate $\hat{\pi}_n$. A similar prior augmentation was studied by Ray and Szabó (2019) for continuous outcomes, and they evaluated the performance of their method empirically when $\hat{\pi}_n$ is trained on the same data as the posterior. Using a point estimate in this fashion is similar to classical correction approaches like targeted maximum likelihood estimation (TMLE) (van der Laan & Rose, 2011; van der Laan & Rubin, 2006). As we will discuss further in Section 6.3, we have found that it is beneficial to instead incorporate the Bayesian uncertainty in the nuisance parameters (in this case the propensity score) by propagating the variation from the initial posterior to the one-step posterior.

# 5 The average treatment effect on the treated

Consider the problem of inferring the causal effect of a binary treatment A on an outcome Y using observational data. Our data takes the form $Z = (X, A, Y)$, where X is a vector of pretreatment covariates. Causal frameworks often posit the existence of counterfactual variables and define causal effects as contrasts between functionals of counterfactual distributions (e.g. Hernán & Robins, 2020; Pearl, 2009). In this case, we can define the pair of variables $(Y^1, Y^0)$ to be the potential outcomes associated with treatment values 1 (treated) and 0 (controls), respectively. Under the consistency assumption (Hernán & Robins, 2020)

$$Y = AY^1 + (1 - A)Y^0,$$

we observe exactly one of $(Y^1, Y^0)$ for each unit in the dataset.

Our target estimand is the *average treatment effect on the treated* (ATT) (Heckman & Robb, 1985; Rubin, 1977)

$$\chi(P) = E_P[Y^1 - Y^0 | A = 1],$$

which quantifies the average causal effect of withholding treatment on the population currently being treated. This can also be formulated without counterfactuals (e.g. Dawid, 2021; Geneletti & Dawid, 2011).

An attractive aspect of the ATT is that it is identified under assumptions that are strictly weaker than those commonly specified for the more well-known average treatment effect $E_P[Y^1 - Y^0]$. More specifically, the identification formula

$$\chi(P) = E_P[Y - E_P[Y | A = 0, X] | A = 1] \qquad (21)$$

can be obtained under the assumptions of *unconfoundedness for controls*

$$Y^0 \perp A | X$$

and *weak overlap*:

$$P(A = 1 | X) < 1$$
with P-probability 1 (Imbens, 2004). Under violations of consistency or unconfoundedness for controls, our results still apply for estimation of the functional defined by (21) provided that weak overlap holds.

Similar to Section 4, it is natural to parameterize P in terms of the propensity score $\pi(x) = P(A = 1 | X = x)$, the conditional outcome density $p_{Y|A,X}(y | a, x)$, and the marginal covariate density $q$. For a single observation $z = (x, a, y)$, the density

$$p(z) = \pi(x)^a (1 - \pi(x))^{1-a} p_{Y|A,X}(y | a, x) q(x)$$

factorizes into the three parameters, so it is particularly convenient and intuitive to specify independent priors.

The functional in (21) has efficient influence function

$$\dot{\chi}_P = \frac{A - \pi(X)}{P[A][1 - \pi(X)]} \{Y - \mu^{(0)}(X)\} - \frac{A}{P[A]} \chi(P), \quad (22)$$

where $\mu^{(0)}(x) = E_P[Y | A = 0, X = x]$ is the outcome regression function for the controls (Chapter 8 of van der Laan & Rose, 2011). We highlight the presence of the marginal probability of treatment assignment $P[A]$, which we can write as $P[A] = Q[\pi(X)]$ with respect to the natural parameterization in terms of $\pi$ and the marginal covariate distribution $Q$.

The behaviour of a posterior on $Q[\pi(X)]$ will depend on the interplay between the models for $\pi$ and $Q$. In the construction of our one-step posterior, we have found it to be beneficial theoretically to replace $Q[\pi(X)]$ with $\tilde{P}[A]$, which is partially due to the fact that the posterior for $\tilde{P}[A]$ is guaranteed to contract to $P_0[A]$ at the rate of $n^{-1/2}$ without requiring additional assumptions on $\pi$ and $Q$. With a slight abuse of notation and terminology, we define the one-step corrected parameter to be

$$\tilde{\chi}(P, \tilde{P}) = \tilde{P} \left[ \frac{A - \pi(X)}{\tilde{P}[A][1 - \pi(X)]} \{Y - \mu^{(0)}(X)\} \right]. \quad (23)$$

We clarify that the $\tilde{P}$ in the denominator of the expression is the same as the one that appears outside of the large square brackets. From a computational perspective, this means that the user does not need to draw an additional set of uniform Dirichlet weights for each posterior sample of $P$. Moreover, (23) only depends on $P$ through $\pi$ and $\mu^{(0)}$, so the user can undertake an initial posterior analysis that is conditional on the covariates and avoid computing the posterior for $Q$.

This example helps to illustrate that our framework offers the flexibility to attain more attractive theoretical properties by making modifications to the correction. Similar modifications are used for classical one-step estimation (e.g. Rotnitzky et al., 2021) with the empirical measure and the Bayesian bootstrap once again playing analogous roles. The corrected parameter defined in (23) also admits an alternative motivation by using the efficient influence function in (22) to define a Bayesian-bootstrap-weighted estimating equation

$$\tilde{P} \left[ \frac{A - \pi(X)}{P[A][1 - \pi(X)]} \{Y - \mu^{(0)}(X)\} - \frac{A}{P[A]} \chi \right] = 0,$$

which is solved by $\chi = \tilde{\chi}(P, \tilde{P})$. The relationship between one-step and estimating equation estimators is discussed in Hines et al. (2022) from a classical perspective.

**Assumption 4** Let $\pi_0$ and $\mu_0^{(0)}$ denote the true parameters associated with $P_0$. There exists a sequence of measurable subsets $(H_n)_n$ of $\mathcal{P}$ satisfying $\Pi(P \in H_n | Z_{1:n}) \xrightarrow{P_0} 1$ such that

(a) $(L_2$-convergence of $\pi$ and $\mu^{(0)})$ There exist numbers $\rho_n, \varepsilon_n \to 0$ such that $\sqrt{n}\rho_n\varepsilon_n \to 0$ and

$$\sup_{P \in H_n} \|\pi - \pi_0\|_{P_0} \le \rho_n \quad \text{and} \quad \sup_{P \in H_n} \|\mu^{(0)} - \mu_0^{(0)}\|_{P_0} \le \varepsilon_n.$$
(b) (Uniform bounding) For all sufficiently large $n$, there exist fixed constants $\delta$, $C > 0$ such that for all $P \in H_n$,

$$ \|\pi\|_{\infty}, \|\pi_0\|_{\infty} < 1 - \delta \quad \text{and} \quad P_0(|Y - \mu^{(0)}(X)| < C) = 1. $$

(c) (Donsker class) The sequences of sets $\{\pi: P \in H_n\}$ and $\{\mu^{(0)}: P \in H_n\}$ are both eventually contained in fixed $P_0$-Donsker classes.

**Theorem 6** Under Assumption 4, the one-step posterior for the ATT functional (21) satisfies the semiparametric BvM theorem.

The proof of Theorem 6 can be found in online supplementary material, Section S8.7. Assumption 4 has a similar form to Assumption 2 for the example in Section 4. In particular, the second-order bias in this case possesses a cross-term structure in terms of $\pi$ and $\mu^{(0)}$. We provide an example of verifying Assumption 4 in online supplementary material, Section S5.2.

Given the conspicuous absence of $\mu^{(1)}(x) = E_P[Y | A = 1, X = x]$ from the definitions and assumptions, it is natural to ask whether one should model the parameters $p_{Y|A,X}(y | 1, x)$ and $p_{Y|A,X}(y | 0, x)$ with independent priors, such that the model for $p_{Y|A,X}(y | 1, x)$ drops out of the analysis for $\tilde{\pi}$, and the posterior for $\mu^{(0)}$ is derived solely from the data on the controls. In the heterogeneous causal effects literature, this separation of the regression model into the two treatment arms is called the *T-learner* approach ('two' learners; Künzel et al., 2019).

In general, we would recommend against doing this and instead advocate for specifying a single regression model that treats $A$ as just another covariate, which is called the *S-learner* approach ('single' learner; Künzel et al., 2019). Provided that Assumption 4 is satisfied, there is no difference in terms of asymptotic performance, but we would expect to obtain finite-sample efficiency gains from the pooling of information across the entire dataset.

An advantage of using the T-learner approach, however, would be to prevent potential model misspecification for $p_{Y|A,X}(y | 1, x)$ from contaminating the posterior for $p_{Y|A,X}(y | 0, x)$, thereby reducing the risk of violating Assumption 4(a). There could also be a finite-sample impact if $\mu_0^{(0)}$ and $\mu_0^{(1)}$ are very different (e.g. $\mu_0^{(0)}$ is much smoother than $\mu_0^{(1)}$ due to high heterogeneity in the conditional treatment effect given $X$), because prior regularization may shrink a joint posterior for $p_{Y|A,X}(y | a, x)$ towards a compromise that is detrimental to precision for $\mu_0^{(0)}$. Ultimately, we would encourage the user to think carefully about their modelling specifications, both in terms of incorporating substantive knowledge and the implications concerning the theory we have provided.

# 6 Simulation studies

## 6.1 Estimating the integrated squared density under a Laplace distribution

Suppose that the data $Z_{1:n}$ are generated i.i.d. from a Laplace(0, 1) distribution, which has density

$$ f_0(z) = \frac{1}{2} \exp(-|z|). $$

The Laplace distribution was chosen because it is nondifferentiable at $z = 0$, and we would expect this to induce oversmoothing in the Dirichlet process Gaussian mixture model from Section 3. We specified our (initial) Bayesian model according to (7) with the standard normal distribution for the base measure $G_0$ and the prior parameter values $\{M = 1, a = 1, b = 1\}$.

Our goal is to estimate the integrated squared density as described in Section 3. To evaluate the performance of the one-step posterior correction, we ran 1,000 independent Monte Carlo trials for each of the two sample sizes $n = 1,000$ and $n = 2,000$. In addition to the uncorrected and corrected posteriors (referred to as 'DPMM' and 'DPMM+1step', respectively), we also include results for the posterior that replaces the Bayesian bootstrap parameter $\tilde{P}$ with the empirical measure $\tilde{P}_n$ in the one-step parameter formulation. We refer to this method as 'DPMM+Emp'. For all methods, point estimation was evaluated using the posterior mean, and intervals were constructed using the 95% central credible region.
<table><caption>**Table 1.** Monte Carlo numerical results comparing the integrated squared density posteriors with and without correction</caption><thead><tr><th>Sample size</th><th>Method</th><th>Bias</th><th>MAE</th><th>RMSE</th><th>Cov (%)</th><th>Int. len.</th></tr></thead><tbody><tr><td rowspan="3">n = 1,000</td><td>DPMM</td><td>-0.0098</td><td>0.0108</td><td>0.0185</td><td>64.9</td><td>0.0290</td></tr><tr><td>DPMM+Emp</td><td>-0.0040</td><td>0.0072</td><td>0.0151</td><td>82.8</td><td>0.0304</td></tr><tr><td>DPMM+1step</td><td>-0.0040</td><td>0.0072</td><td>0.0151</td><td>93.7</td><td>0.0434</td></tr><tr><td rowspan="3">n = 2,000</td><td>DPMM</td><td>-0.0068</td><td>0.0073</td><td>0.0138</td><td>67.4</td><td>0.0210</td></tr><tr><td>DPMM+Emp</td><td>-0.0025</td><td>0.0048</td><td>0.0116</td><td>87.5</td><td>0.0234</td></tr><tr><td>DPMM+1step</td><td>-0.0025</td><td>0.0048</td><td>0.0116</td><td>95.8</td><td>0.0324</td></tr></tbody></table>

Note. MAE = median absolute error; RMSE = root mean squared error; Cov = coverage; Int. len. = average 95% central credible interval length.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Comparison of sampling distributions for the Bayesian point estimators of the integrated squared density, illustrated by histogram/density plots: *n* = 1,000 (left), and *n* = 2,000 (right). The true value of the functional is shown by the black vertical dashed line.</figcaption></figure>

The numerical results can be found in Table 1. First, we can see that the one-step posterior improved on the initial posterior across all point estimation metrics. This was achieved through substantial bias reduction, which is illustrated visually by the histogram/density plots of the posterior mean in Figure 2. The point estimation results for ‘DPMM+Emp’ and ‘DPMM+1step’ are the same because the conditional expectation of $\hat{P}$ given $Z_{1:n}$ is equal to $P_n$.

The improvements in coverage are even more striking. For the larger sample size of *n* = 2,000, the one-step posterior was able to attain nominal coverage despite the fact that the initial posterior only achieved 67.4%. The poor performance of the initial posterior can be partially attributed to the aforementioned bias, which is of a similar order of magnitude to the average interval length. However, it also appears that the interval lengths are too short. This could have been caused by an overcurved likelihood or residual prior regularization, or both. In any case, it seems unlikely that the initial posterior satisfies the semiparametric BvM theorem for this Laplace example. As expected, the ‘DPMM+Emp’ posterior also substantially undercovers despite the bias reduction. This further emphasizes our point in Section 2.2 that the Bayesian bootstrap is crucial for correcting the shape of the posterior.

For further context, we have included a ‘before and after’ scatter plot in Figure 3 that compares the samples from both posteriors for a single iteration of the experiment (*n* = 2,000). This plot
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Scatter plot comparing the initial posterior samples with their one-step corrected counterparts. We have also displayed the identity line (—) and the true value of the estimand (---).</figcaption></figure>

illustrates that, despite the randomness of the correction terms, the draws from the initial posterior tend to be pulled in the direction of the true functional value.

## 6.2 Missing data

To evaluate performance for the missing data example of Example 1 and Section 4, we use a data-generating process from Kennedy (2023). The covariates X are uniform on [-1, 1], the propensity score is $\pi(x) = 0.5 + 0.4 \text{sign}(x)$, and the outcome regression function is defined by

$$m(x) = \begin{cases} \frac{(x+2)^2}{2} & \text{if } -1 \le x < -0.5 \\ \frac{x}{2} + 0.875 & \text{if } -0.5 \le x < 0 \\ -5(x-0.2)^2 + 1.075 & \text{if } 0 \le x < 0.5 \\ x + 0.125 & \text{if } 0.5 \le x \le 1. \end{cases}$$

Conditional on $X=x$, the outcomes are generated by $Y \sim N(m(x), (0.2 - 0.1 \cos(2\pi x))^2)$. This piecewise polynomial structure of $m(x)$ originated from Györfi et al. (2002), and it was selected by Kennedy (2023) to illustrate the risks of under/oversmoothing faced by nonparametric regression methods under missingness.

For our initial posterior, we modelled both the conditional outcome distribution and propensity score with BART (Chipman et al., 2010). We used only the default hyperparameters in the R BART package implementations with no cross-validation. The propensity score model used the pbart procedure, which specifies a probit link. We drew 5,000 posterior samples for each iteration, after burning 1,000 samples for initiation. To perform standard Bayesian inference for the functional, we also specify the Bayesian bootstrap for the marginal covariate distribution, which has become a
<table><caption>**Table 2.** Monte Carlo numerical results comparing the outcome mean estimation procedures</caption><thead><tr><th>Sample size</th><th>Method</th><th>Bias</th><th>MAE</th><th>RMSE</th><th>Cov (%)</th><th>Int. len.</th></tr></thead><tbody><tr><td rowspan="5">n = 500</td><td>Standard BART</td><td>0.0062</td><td>0.0182</td><td>0.0275</td><td>91.0</td><td>0.0933</td></tr><tr><td>Plug-in+EIF</td><td>0.0034</td><td>0.0184</td><td>0.0273</td><td>84.0</td><td>0.0773</td></tr><tr><td>Plug-in+BB</td><td>0.0034</td><td>0.0184</td><td>0.0273</td><td>83.5</td><td>0.0781</td></tr><tr><td>Antonelli</td><td>0.0026</td><td>0.0186</td><td>0.0273</td><td>92.1</td><td>0.0990</td></tr><tr><td>One-step posterior</td><td>0.0026</td><td>0.0186</td><td>0.0273</td><td>94.7</td><td>0.1118</td></tr><tr><td rowspan="5">n = 1,000</td><td>Standard BART</td><td>0.0033</td><td>0.0113</td><td>0.0177</td><td>93.5</td><td>0.0666</td></tr><tr><td>Plug-in+EIF</td><td>0.0017</td><td>0.0113</td><td>0.0177</td><td>88.1</td><td>0.0568</td></tr><tr><td>Plug-in+BB</td><td>0.0018</td><td>0.0112</td><td>0.0177</td><td>88.4</td><td>0.0572</td></tr><tr><td>Antonelli</td><td>0.0014</td><td>0.0112</td><td>0.0177</td><td>94.1</td><td>0.0684</td></tr><tr><td>One-step posterior</td><td>0.0014</td><td>0.0112</td><td>0.0177</td><td>96.1</td><td>0.0761</td></tr><tr><td rowspan="5">n = 2,000</td><td>Standard BART</td><td>0.0016</td><td>0.0085</td><td>0.0124</td><td>94.3</td><td>0.0470</td></tr><tr><td>Plug-in+EIF</td><td>0.0008</td><td>0.0085</td><td>0.0124</td><td>90.4</td><td>0.0418</td></tr><tr><td>Plug-in+BB</td><td>0.0008</td><td>0.0085</td><td>0.0124</td><td>90.3</td><td>0.0420</td></tr><tr><td>Antonelli</td><td>0.0006</td><td>0.0085</td><td>0.0124</td><td>94.5</td><td>0.0480</td></tr><tr><td>One-step posterior</td><td>0.0006</td><td>0.0085</td><td>0.0124</td><td>96.5</td><td>0.0520</td></tr></tbody></table>

Note. MAE = median absolute error; RMSE = root mean squared error; Cov = coverage; Int. len. = average 95% interval length.

popular default choice for Bayesian causal inference (Linero & Antonelli, 2022; Ray & van der Vaart, 2020).

Our first comparator is classical one-step estimation. In each iteration, we computed the posterior means ($\hat{m}$, $\hat{\pi}$) from the BART posteriors and then plugged them into the one-step estimator. For interval estimation, we used the sample variance of the estimated efficient influence function—as described in Section 2.3—combined with a Gaussian approximation to the sampling distribution. We call this ‘Plug-in+EIF’. Similarly, we used the posterior means ($\hat{m}$, $\hat{\pi}$) to implement the Dirac mass one-step posterior idea from Section 2.3, which we call ‘Plug-in+BB’. The last comparator is the proposal from Antonelli et al. (2022) with the same BART posteriors (see online supplementary material, Section S4 for more detailed descriptions). This method requires computing B bootstrap samples for every posterior sample. We closely followed the implementation from https://github.com/jantonelli111/DoublyRobustHD, where B was set to 100. In accordance with this implementation, the samples of the propensity score were trimmed to lie in [0.05, 0.95].

Our numerical results can be found in Table 2. We ran 2,000 Monte Carlo trials for each of the sample sizes $n = 500, 1,000, 2,000$. For the Bayesian methods, we evaluated point estimation with the posterior mean and interval estimation with the central 95% credible intervals. The results for point estimation are very similar across the board, so we focus our attention on interval estimation. First, we see that the benchmark ‘Standard BART’ undercovers at all sample sizes despite the fact that the magnitude of bias is small relative to the interval lengths in this case. The coverage does improve as the sample size increases however. More strikingly, the two plug-in methods perform very poorly, even coming behind ‘Standard BART’. This strongly suggests that the ignored second-order remainder terms retain a significant impact, and the interval lengths are too short as a result. As expected, there is very little difference between the two plug-in methods.

The ‘Antonelli’ method improves on both ‘Standard BART’ and the plug-in approaches, but there is still slight undercoverage. On the other hand, the one-step posterior was able to obtain nominal coverage for $n = 1,000$ and $n = 2,000$. In online supplementary material, Section S4, we show that the variance estimator of ‘Antonelli’ is always bounded above by the one-step posterior variance, which is confirmed by our numerical results here. Compared to ‘Antonelli’, the one-step posterior interval lengths were 12.9% longer on average at $n = 500$, then 11.3% longer at $n = 1,000$, and 8.3% longer at $n = 2,000$. This suggests, reassuringly, that there can be
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Breakdown of procedure run times averaged across 10 iterations each. The times for the Antonelli variance estimation and one-step posterior sampling do not include the initial BART posterior sampling.</figcaption></figure>

significant differences between the two methods for smaller sample sizes, but they converge together in the asymptotic regime.

Finally, we provide a breakdown of computational demands in Figure 4. As benchmarks, we have displayed the run times for the outcome regression and propensity score BART models. Then we have the times for the variance estimation procedure of Antonelli et al. (2022), which do not include the BART posterior sampling. We see that the Antonelli et al. (2022) method is substantially more computationally intensive than one-step posterior sampling. This is because they must compute B bootstrap samples for every posterior sample, while we only bootstrap one-to-one. More details can be found in online supplementary material, Section S4.

## 6.3 The ACIC 2016 competition

The causal data analysis competition ‘Is Your SATT Where It’s At?’ was held at the 2016 Atlantic Causal Inference Conference with the goal of providing a comprehensive, objective comparison of the available methods for estimating causal effects. A total of 24 submissions were compared in the ‘black box’ section of the competition, which assessed the automated performance of the methods across 77 different simulation settings. Detailed descriptions of the competition design, motivations and findings can be found in Dorie et al. (2019).

Each dataset $Z_{1:n} = (Z_1, ..., Z_n)$ had a sample size of $n = 4, 802$ with the same data structure as Section 5; that is, $Z_i = (X_i, A_i, Y_i)$ for $i = 1, ..., n$. The set of covariates $(X_1, ..., X_n)$ was fixed across datasets—they were taken from a real study that examined the impact of birth weight on a child’s IQ and contained a mixture of categorical, binary, continuous and count variables (58 variables in total). Conditional on the covariates, the potential outcomes $(Y_i^1, Y_i^0)$ and treatment assignments were generated i.i.d. according to the configuration of simulation parameters, which included the degree of nonlinearity of the regression functions, the treated percentage, and the treatment effect heterogeneity. One hundred independent datasets were generated from each of the 77 settings, resulting in a total of 7,700 datasets.

The target estimand is the *sample average treatment effect on the treated* (SATT)

$$ \beta = \frac{\mathbb{P}_n[A(Y^1 - Y^0)]}{\mathbb{P}_n[A]}, \qquad (24) $$

which quantifies the average treatment effect among those treated in the dataset. But the SATT cannot be directly inferred in a semiparametric framework because we never observe $AY^0$ when
A = 1. An obvious alternative is to use the estimates for the (population) ATT as described in Section 5, but this would likely produce substantially conservative intervals for the SATT (Imbens, 2004). As a partial remedy, we will also consider estimates derived from the one-step posterior targeting a new estimand that we call the *average covariate-conditional treatment effect on the treated* (ACTT)

$$ \theta(P) = \frac{\mathbb{Q}_n[\pi(X)\{\mu^{(1)}(X) - \mu^{(0)}(X)\}]}{\mathbb{Q}_n[\pi(X)]}, $$

which replaces the unknown population covariate distribution $Q$ in the ATT identification formula (21) with the empirical covariate distribution

$$ \mathbb{Q}_n = \frac{1}{n} \sum_{i=1}^{n} \delta_{X_i}. $$

We refer the reader to online supplementary material, Section S6 for further information about the ACTT, including details on implementing the one-step posterior correction.

Like the analysis in Section 6.2, we used BART with the default hyperparameters from the R BART package to model the conditional outcome distribution and propensity score, and the Bayesian bootstrap to model the marginal covariate distribution. As discussed at the end of Section 5, we fitted the outcome regression model on both treatment arms together to pool information. For each dataset, we obtained 4,000 posterior samples aggregated from 4 independent chains after discarding 2,000 burn-in samples. The posterior mean and the central 95% credible interval were used for point and interval estimation, respectively.

We highlight two top performers from the competition (descriptions of the remaining submissions can be found on p. 62 of Dorie et al., 2019). The ‘BART on Pscore’ method (Hahn et al., 2020) fitted a BART outcome regression model that included an estimate of the propensity score (from a cross-validated probit BART model) as a splitting variable. A posterior for $\beta$ was then derived by replacing $(Y_i^1, Y_i^0)$ with $(\mu^{(1)}(X_i), \mu^{(0)}(X_i))$ in (24) (Hill, 2011). The argument put forth by Hahn et al. (2020) was that their augmentation helps to mitigate regularization bias when estimating the outcome regression function because the propensity score is likely to be correlated with $\mu_0^{(0)}$ in practice, e.g. a physician may assign treatment more readily to patients with worse risk factors. Accordingly, the propensity score should be an informative transformation of the covariates for predicting the outcome, which makes it easier for the BART model to approximate the truth with finite samples.

The ‘BART+TMLE’ method employed a correction procedure similar to that of the one-step estimator. TMLE (van der Laan & Rose, 2011; van der Laan & Rubin, 2006) makes adjustments in the model space rather than the parameter space; an initial estimate $\tilde{P}$ is mapped to a targeted estimate $\tilde{P}$ such that $\mathbb{P}_n[\tilde{\chi}_P] \approx 0$ for the population ATT. Consequently, the TMLE plug-in $\chi(\tilde{P})$ is approximately equal to its own one-step estimator and its variance can be estimated by $\frac{1}{n} \mathbb{P}_n[(\tilde{\chi}_P - \mathbb{P}_n[\tilde{\chi}_P])^2]$ provided that $\tilde{P}$ satisfies the requisite conditions for asymptotic efficiency (see online supplementary material, Section S1.2). In this case, the initial estimates of the outcome regression function and the propensity score were both derived from BART.

Figure 5 compares the standardized point estimation errors across all settings. The top performer from the competition submissions was ‘BART on Pscore’. In particular, it improved on ‘BART MChains’, which was the same procedure aside from the propensity score augmentation (‘MChains’ stands for multiple chains). We emphasize again, however, that the ‘BART on Pscore’ augmentation does not target any particular low-dimensional functional, and the Bayesian update will instead use the estimated propensity score to help find a favourable bias-variance trade-off for the whole outcome regression function. In contrast, our one-step posteriors used the propensity score model to debias specifically in the direction of the ATT, which is perhaps why they were able to present further improvements in reducing the point estimation errors. However, as Hahn et al. (2020) pointed out, there is the option of combining ‘BART on Pscore’ with a targeted approach like TMLE, and by extension, we can also combine it with the one-step posterior correction to reap the benefits of both approaches. For additional context,
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparison of standardized point estimation errors for ACIC 2016. The competition submissions are shown in black. Root mean squared errors are shown by the blue dots; the median, lower quartile and upper quartile absolute errors across the 7,700 datasets are shown in orange.</figcaption></figure>

we have included the results for the BART ACTT posterior without the one-step correction [BART plug-in (ACTT)] and it is clear that the one-step ACTT posterior performs substantially better.

The results for interval estimation can be found in Figure 6. None of the competition submissions were able to attain the nominal coverage level of 95% across all 7,700 datasets, with ‘BART+TMLE’ coming the closest. We recall that the intervals for ‘BART+TMLE’ were constructed for the population ATT rather than the SATT, so we would expect them to be conservative. As the direct comparator, the one-step ATT posterior was able to attain a conservative level of coverage, albeit with slightly wider intervals. As expected, the one-step ACTT posterior had shorter interval lengths than both methods, but it was still able to attain nominal coverage. Taking both coverage and interval lengths into account, the one-step ACTT posterior was the clear winner in this category. As before, we have included the results for the uncorrected ACTT posterior, which performed slightly worse than the other plug-in BART methods like ‘BART on Pscore’ and ‘BART Mchains’.

Earlier, we described how the variance of the TMLE estimator could be estimated by $n^{-1}\mathbb{P}_n[(\hat{\chi}_P - \mathbb{P}_n[\hat{\chi}_P])^2]$, which is the sample variance of the estimated efficient influence function divided by the sample size. As discussed in Section 2.3, other influence-function-based methods like one-step estimation take the same approach but with $\hat{P}$ replaced by the original estimate $\hat{P}$. This variance estimate is similarly based on a first-order asymptotic approximation that assumes the second-order remainder terms to be negligible. Like the results in Section 6.2, we see that the intervals derived from these variance estimates tend to undercover in practice.

In conclusion, the one-step posterior is a much more natural and effective way to leverage Bayesian procedures for semiparametric inference than classical approaches like one-step estimation and TMLE. We expect there to be the potential for further performance gains for one-step BART posteriors if the user employs cross-validation to select hyperparameters or augments the outcome regression model with an estimated propensity score augmentation like ‘BART on Pscore’.

# 7 Discussion

The holistic nature of Bayesian inference—automatic uncertainty quantification, marginalizing nuisance parameters, inference for all estimands in one stroke—is often cited as one of its key
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Comparison of coverage and average standardized interval lengths for ACIC 2016. The competition submissions are shown in black. The dots are averages across all 7,700 datasets, while the bars show the upper and lower quartiles for the 77 settings. The dashed horizontal line indicates the 95% coverage level.</figcaption></figure>

strengths. This wisdom is subverted in the semiparametric setting. A Bayesian user has the unenviable choice of either tailoring their model towards a particular estimand at the expense of all others or run the risk of obtaining biased, uncalibrated inference. Our one-step posterior correction allows us to circumvent this problem by targeting estimands after the initial model has been fitted.

A potential obstacle for our approach is the requirement that the efficient influence function can be computed across the whole model. This could be particularly challenging for separated semi-parametric models, because the efficient influence function is defined in terms of orthogonal projections (van der Vaart, 1998) and may not admit a tractable expression. However, with recent developments in numerical procedures for approximating efficient influence functions (Carone et al., 2019; Jordan et al., 2022), we anticipate that it will be possible to fully automate the one-step posterior correction in future work.

Finally, we discuss a utility that might be particularly appealing to practicing Bayesians. Typically, a nonparametric Bayesian model is carefully constructed based on a number of considerations, such as computational scalability, prior information, interpretability etc. Our framework then provides a form of sensitivity analysis, in which the initial posterior could be compared with the one-step posterior. If the two posteriors are in broad agreement (e.g. both posteriors indicate a significant causal effect in the same direction), then this could be interpreted as a reinforcement of the confidence in the analysis. On the other hand, if there is a disparity between the two posteriors, then this indicates potential limitations in the initial posterior, such as an overly influential prior that contradicts the evidence in the data. In this case, the user might prefer to report only the one-step posterior analysis because of its small bias property.

## Acknowledgments

E.F. was an employee of Novo Nordisk during part of this work. The authors thank Edward Kennedy, Kolyan Ray, the joint editor, the associate editor, and two anonymous reviewers for their helpful input that have improved the manuscript.
*Conflicts of interest:* None declared.

**Funding**

A.Y. received funding from Novo Nordisk. C.H. was supported by the Alan Turing Institute, the Li Ka Shing Foundation, the Medical Research Council, the EPSRC through the Bayes4Health grant EP/R018561/1, AI for Science and Government UKRI, and the U.K. Engineering and Physical Sciences Research Council. J.R. has received funding from the European Research Council (ERC) under the Horizon 2020 research and innovation programme (grant agreement no. 834175).

**Data availability**

All datasets used in this paper are publicly available. Code for reproducing the simulation results in the paper can be found at https://github.com/andrew-yiu/Semiparametric-posterior-corrections.

**Supplementary material**

Supplementary material is available online at *Journal of the Royal Statistical Society: Series B*.

**References**

Alquier P., & Ridgway J. (2019). Concentration of tempered posteriors and of their variational approximations. *Annals of Statistics*, 48(3), 1475-1497. https://doi.org/10.1214/19-AOS1855

Antonelli J., Papadogeorgou G., & Dominici F. (2022). Causal inference in high dimensions: A marriage between Bayesian modeling and good frequentist properties. *Biometrics*, 78(1), 100-114. https://doi.org/10.1111/biom.v78.1

Berger J. (1985). *Statistical decision theory and Bayesian analysis*. Springer-Verlag.

Berger J. O., Liseo B., & Wolpert R. L. (1999). Integrated likelihood methods for eliminating nuisance parameters. *Statistical Science*, 14(1), 1-28. https://doi.org/10.1214/ss/1009211804

Bickel P., & Ritov Y. (1988). Estimating integrated squared density derivatives: Sharp best order of convergence estimates. *Sankhya*, 50, 381-393.

Bickel P., & Ritov Y. (2003). Nonparametric estimators which can be “plugged-in”?. *Annals of Statistics*, 31(4), 1033-1053. https://doi.org/10.1214/ aos/1059655904

Birgé L., & Massart P. (1995). Estimation of integral functionals of a density. *Annals of Statistics*, 23, 11-29. https://doi.org/10.1214/ aos/1176324452

Breiman L. (2001). Random forests. *Machine Learning*, 45(1), 5-32. https://doi.org/10.1023/A:1010933404324

Buja A., Brown L., Berk R., George E., Pitkin E., Traskin M., Zhang K., & Zhao L. (2019a). Models as approximations I: Consequences illustrated with linear regression. *Statistical Science*, 34, 523-544. https://doi.org/10.1214/18-STS693

Buja A., Brown L., Kuchibhotla A., Berk R., George E., & Zhao L. (2019b). Models as approximations II: A model-free theory of parametric regression. *Statistical Science*, 34, 545-565. https://doi.org/10.1214/18-STS694

Carone M., Luedtke A., & van der Laan M. (2019). Toward computerized efficient estimation in infinite-dimensional models. *Journal of the American Statistical Association*, 114(527), 1174-1190. https://doi.org/10.1080/01621459.2018.1482752

Castillo I. (2012). A semiparametric Bernstein-von Mises theorem for Gaussian process priors. *Probability Theory and Related Fields*, 152(1-2), 53-99. https://doi.org/10.1007/s00440-010-0316-5

Castillo I., & Nickl R. (2014). On the Bernstein-von Mises phenomenon for nonparametric Bayes procedures. *Annals of Statistics*, 42(5), 1941-1969. https://doi.org/10.1214/14-AOS1246.

Castillo I., & Rousseau J. (2015). A Bernstein-von Mises theorem for smooth functionals in semiparametric models. *Annals of Statistics*, 43, 2353-2383. https://doi.org/10.1214/15-AOS1336

Cheng G., & Kosorok M. (2008a). General frequentist properties of the posterior profile distribution. *Annals of Statistics*, 36, 1819-1853. https://doi.org/10.1214/07-AOS536

Cheng G., & Kosorok M. (2008b). Higher order semiparametric frequentist inference with the profile sampler. *Annals of Statistics*, 36, 1786-1818. https://doi.org/10.1214/07-AOS523

Chernozhukov V., Chetverikov D., Demirer M., Duflo E., Hansen C., Newey W., & Robins J. (2018). Double/debiased machine learning for treatment and structural parameters. *The Econometrics Journal*, 21(1), C1-C68. https://doi.org/10.1111/ectj.12097

Chipman H., George E., & McCulloch R. (2007). Bayesian ensemble learning. *Neural Information Processing Systems*, 19, 265-272. https://doi.org/10.7551/mitpress/7503.003.0038
Chipman H., George E., & McCulloch R. (2010). BART: Bayesian additive regression trees. *Annals of Applied Statistics*, 4(1), 266–298. https://doi.org/10.1214/09-AOAS285

Dawid A. (2021). Decision-theoretic foundations for statistical causality. *Journal of Causal Inference*, 9(1), 39–77. https://doi.org/10.1515/jci-2020-0008

Diaconis P., & Freedman D. (1986). On the consistency of Bayes estimates. *Annals of Statistics*, 14(1), 1–26. https://doi.org/10.1214/aos/1176349830.

Donnet S., Rivoirard V., Rousseau J., & Scricciolo C. (2018). Posterior concentration rates for empirical Bayes procedures with applications to Dirichlet process mixtures. *Bernoulli*, 24(1), 231–256. https://doi.org/10.3150/16-BEJ872

Dorie V., Hill J., Shalit U., Scott M., & Cervone D. (2019). Automated versus do-it-yourself methods for causal inference: Lessons learned from a data analysis competition. *Statistical Science*, 34, 43–68. https://doi.org/10.1214/18-STS667

Efron B., & Hastie T. (2016). *Computer age statistical inference*. Cambridge University Press.

Ferguson T. (1973). A Bayesian analysis of some nonparametric problems. *Annals of Statistics*, 2, 209–230. https://doi.org/10.1214/aos/1176342360

Fisher A., & Kennedy E. (2021). Visually communicating and teaching intuition for influence functions. *The American Statistician*, 70(2), 162–172. https://doi.org/10.1080/00031305.2020.1717620

Geneletti S., & Dawid A. (2011). Defining and identifying the effect of treatment on the treated. In P. Illari, F. Russo, & J. Williamson (Eds.), *Causality in the sciences* (pp. 728–749). Oxford University Press.

Ghosal S., & van der Vaart A. (2007). Convergence rates for noniid observations. *Annals of Statistics*, 35, 192–223. https://doi.org/10.1214/aos/1016218228

Ghosal S., & van der Vaart A. (2017). *Fundamentals of nonparametric Bayesian inference*. Cambridge University Press.

Ghosh J., & Ramamoorthi R. (2003). *Bayesian nonparametrics*. Springer.

Giné E., & Nickl R. (2008). A simple adaptive estimator of the integrated square of a density. *Bernoulli*, 14, 47–61. https://doi.org/10.3150/07-BEJ110

Giné E., & Nickl R. (2016). *Mathematical foundations of infinite-dimensional statistical models*. Cambridge University Press.

Gruber S., & van der Laan M. (2019). Comment on “Automated versus do-it-yourself methods for causal inference: Lessons learned from a data analysis competition”?. *Statistical Science*, 34(1), 82–85. https://doi.org/10.1214/18-STS689

Györfi L., Kohler M., Krzyzak A., & Walk M. (2002). *A distribution-free theory of nonparametric regression*. Springer.

Hahn P., Murray J., & Carvalho C. (2020). Bayesian regression tree models for causal inference: Regularization, confounding, and heterogeneous effects. *Bayesian Analysis*, 15(3), 965–1056. https://doi.org/10.1214/19-BA1195

Heckman J., & Robb R. (1985). Alternative methods for evaluating the impact of interventions: An overview. *Journal of Econometrics*, 30(1–2), 239–267. https://doi.org/10.1016/0304-4076(85)90139-3

Hernán M. A., & Robins J. (2020). *Causal inference: What if*. Boca Raton: Chapman & Hall/CRC.

Hill J. (2011). Bayesian nonparametric modeling for causal inference. *Journal of Computational and Graphical Statistics*, 20(1), 217–240. https://doi.org/10.1198/jcgs.2010.08162

Hines O., Dukes O., Diaz-Ordaz K., & Vansteelandt S. (2022). Demystifying statistical learning based on efficient influence functions. *The American Statistician*, 76(3), 292–304. https://doi.org/10.1080/00031305.2021.2021984

Imbens G. (2004). Nonparametric estimation of average treatment effects under exogeneity: A review. *Review of Economics and Statistics*, 86(1), 4–29. https://doi.org/10.1162/003465304323023651

Jordan M., Ghahramani Z., Jaakkola Z., & Saul L. (1999). An introduction to variational methods for graphical models. *Machine Learning*, 37(2), 183–233. https://doi.org/10.1023/A:1007665907178

Jordan M., Wang Y., & Zhou A. (2022). Empirical Gateaux derivative for causal inference. *Advances in Neural Information Processing Systems*, 35, 8512–8525. https://proceedings.neurips.cc/paper_files/paper/2022/hash/3848fef259495bfd04d60cdc5c1b4db7-Abstract-Conference.html

Kang J., & Schafer J. (2007). Demystifying double robustness: A comparison of alternative strategies for estimating a population mean from incomplete data (with discussion). *Statistical Science*, 22, 523–539. https://doi.org/10.1214/07-STS227

Kennedy E. (2023). Towards optimal doubly robust estimation of heterogeneous causal effects. *Electronic Journal of Statistics*, 17(2), 3008–3049. https://doi.org/10.1214/23-EJS2157

Kennedy E., Bonvini M., & Mishler A. (2021). Comment on “statistical modeling: The two cultures”? by Leo Breiman. *Observational Studies*, 7(1), 145–156. https://doi.org/10.1353/obs.2021.0001

Kruijer W., Rousseau J., & van der Vaart A. (2010). Adaptive Bayesian density estimation with location-scale mixtures. *Electronic Journal of Statistics*, 4, 1225–1257. https://doi.org/10.1214/10-EJS584
Künzel S., Sekhon J., Bickel P., & Yu B. (2019). Metalearners for estimating heterogeneous treatment effects using machine learning. *Proceedings of the National Academy of Sciences*, 116(10), 4156-4165. https://doi.org/10.1073/pnas.1804597116

Laurent B. (1996). Efficient estimation of integral functionals of a density. *Annals of Statistics*, 24(2), 659-681. https://doi.org/10.1214/aoa/1032894458

Le Cam L. (1986). *Asymptotic methods in statistical decision theory*. Springer-Verlag.

Linero A., & Antonelli J. (2022). The how and why of Bayesian nonparametric causal inference. *Wiley Interdisciplinary Reviews: Computational Statistics*, 15, e1583. https://doi.org/10.1002/wics.1583

Little R. (2011). Calibrated Bayes, for statistics in general, and missing data in particular. *Statistical Science*, 26, 162-174. https://doi.org/10.1214/10-STS318

Maugis-Rabusseau C., & Michel B. (2013). Adaptive density estimation using finite Gaussian mixtures. *ESAIM: P&S*, 17, 698-724. https://doi.org/10.1051/ps/2012018

Morgan S., & Winship C. (2007). *Counterfactuals and causal inference: Methods and principles for social research*. Cambridge University Press.

Murphy S., & van der Vaart A. (2000). On profile likelihood. *Journal of the American Statistical Association*, 95(450), 449-465. https://doi.org/10.1080/01621459.2000.10474219

Newey W., Hsieh F., & Robins J. (1998). *Undersmoothing and bias corrected functional estimation*. MIT Working Paper.

Ogburn E., & Shpitser I. (2021). Causal modelling: The two cultures. *Observational Studies*, 7(1), 179-183. https://doi.org/10.1353/obs.2021.0006

Pearl J. (2009). Causal inference in statistics: An overview. *Statistics Surveys*, 3, 96-146. https://doi.org/10.1214/09-SS057

Pfanzagl J. (1982). *Contributions to a general asymptotic statistical theory*. Springer.

Prakasa Rao B. (1983). *Nonparametric functional estimation*. Academic Press.

Ray K., & Szabó B. (2019). Debiased Bayesian inference for average treatment effects. *Advances in Neural Information Processing Systems*, 33, 11929-11939. https://proceedings.neurips.cc/paper_files/paper/2019/hash/342285bb2a8cadef22f667eeb6a63732-Abstract.html

Ray K., & van der Vaart A. (2020). Semiparametric Bayesian causal inference. *Annals of Statistics*, 48(5), 2999-3020. https://doi.org/10.1214/19-AOS1919

Ray K., & van der Vaart A. (2021). On the Bernstein-von Mises theorem for the Dirichlet process. *Electronic Journal of Statistics*, 15(1), 2224-2246. https://doi.org/10.1214/21-EJS1821

Rivoirard V., & Rousseau J. (2012). Bernstein-von Mises theorem for linear functionals of the density. *Annals of Statistics*, 40(3), 1489-1523. https://doi.org/10.1214/12-AOS1004

Robins J., Hernán M., & Wasserman L. (2015). Discussion of: On Bayesian estimation of marginal structural models. *Biometrics*, 71(2), 296-299. https://doi.org/10.1111/biom.12273

Robins J., Li L., Tchetgen Tchetgen E., & van der Vaart A. (2017). Minimax estimation of a functional on a structured high-dimensional model. *Annals of Statistics*, 45(5), 1951-1987. https://doi.org/10.1214/16-AOS1515

Robins J., & Ritov Y. (1997). Toward a curse of dimensionality appropriate (CODA) asymptotic theory for semi-parametric models. *Statistics in Medicine*, 16, 285-319. https://doi.org/10.1002/(SICI)1097-0258(19970215)16:3<1.0.CO;2-T

Robins J., & Rotnitzky A. (1992). Recovery of information and adjustment for dependent censoring using surrogate markers. In *AIDS epidemiology-methodological issues* (pp. 297-331). Birkhäuser.

Rotnitzky A., Smucler E., & Robins J. (2021). Characterization of parameters with a mixed bias property. *Biometrika*, 108(1), 231-238. https://doi.org/10.1093/biomet/asaa054

Rousseau J. (2016). On the frequentist properties of Bayesian nonparametric methods. *Annual Review of Statistics and Its Application*, 3(1), 211-231. https://doi.org/10.1146/statistics.2016.3 issue-1

Rubin D. (1976). Inference and missing data. *Biometrika*, 63(3), 581-592. https://doi.org/10.1093/biomet/63.3.581

Rubin D. (1977). Assignment to treatment group on the basis of a covariate. *Journal of Educational and Behavioral Statistics: A Quarterly Publication Sponsored by the American Educational Research Association and the American Statistical Association*, 2(1), 1-26. https://doi.org/10.3102/10769986002001001

Rubin D. (1981). The Bayesian bootstrap. *Annals of Statistics*, 9(1), 130-134. https://doi.org/10.1214/aoa/1176345338

Rubin D. (1984). Bayesianly justifiable and relevant frequency calculations for the applied statistician. *Annals of Statistics*, 12(4), 1151-1172. https://doi.org/10.1214/aoa/1176346785

Scricciolo C. (2014). Adaptive Bayesian density estimation in $L_p$-metrics with Pitman-Yor or normalized inverse-Gaussian process kernel mixtures. *Bayesian Analysis*, 9(2), 475-520. https://doi.org/10.1214/14-BA863

Shen W., Tokdar S., & Ghosal S. (2013). Adaptive Bayesian multivariate density estimation with Dirichlet mixtures. *Biometrika*, 100(3), 623-640. https://doi.org/10.1093/biomet/ast015
Sulem S., Rivoirard V., & Rousseau J. (2023). Scalable and adaptive variational Bayes methods for Hawkes processes. arXiv, arXiv:2212.00293, preprint: not peer reviewed.

van der Laan M., & Robins J. (2003). *Unified methods for censored longitudinal data and causality*. Springer-Verlag.

van der Laan M., & Rose S. (2011). *Targeted learning*. Springer-Verlag.

van der Laan M., & Rose S. (2018). *Targeted learning in data science*. Springer-Verlag.

van der Laan M., & Rubin D. (2006). Targeted maximum likelihood learning. *International Journal of Biostatistics*, 2(1). https://doi.org/10.2202/1557-4679.1043

van der Laan M., Wang Z., & van der Laan L. (2021). 'Higher Order Targeted Maximum Likelihood Estimation', arXiv, arXiv:2101.06290, https://doi.org/10.48550/arXiv.2101.06290, preprint: not peer reviewed.

van der Vaart A. (1991). On differentiable functionals. *Annals of Statistics*, 19, 178–204. https://doi.org/10.1214/ aos/1176347976

van der Vaart A. (1992). Asymptotic linearity of minimax estimators. *Statistica Neerlandica*, 46(2–3), 179–194. https://doi.org/10.1111/stan.1992.46 issue-2-3

van der Vaart A. (1998). *Asymptotic statistics*. Cambridge University Press.

van der Vaart A. (2014). Higher order tangent spaces and influence functions. *Statistical Science*, 29(4), 679–686. https://doi.org/10.1214/14-STS478

van der Vaart A., & Wellner J. (1996). *Weak convergence and empirical processes*. Springer.

van der Vaart A., & Wellner J. (2023). *Weak convergence and empirical processes* (2nd ed.). Springer.

Vansteelandt S. (2021). Statistical modelling in the age of data science. *Observational Studies*, 7(1), 217–228. https://doi.org/10.1353/obs.2021.0013

Vansteelandt S., & Dukes O. (2022). Assumption-lean inference for generalised linear model parameters (with discussion). *Journal of the Royal Statistical Society, Series B*, 84(3), 657–685. https://doi.org/10.1111/rssb. 12504

Zhang F., & Gao C. (2020). Convergence rates of variational posterior distributions. *Annals of Statistics*, 48(4), 2180–2207. https://doi.org/10.1214/19-AOS1883