Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
58,702 characters · 9 sections · 77 citation commands
Robust Semiparametric Inference for Bayesian Additive Regression Trees
The Bayesian additive regression trees (BART) is one of the most powerful statistical methods. Its empirical success in predictive modeling has been demonstrated in numerous studies that uncover complex nonlinear relationships in high-dimensional data. Additionally, in the context of nonparametric function estimation, the BART achieves the adaptive near rate-minimax posterior contraction, as shown by linero2018bart,rockova2020posterior,rockova2020bart,jeong2023art. Beyond prediction, BART has also shown strong performance in semiparametric inference for causal parameters of interest dorie2019auto,hahn2020bart. However, a theoretical justification for semiparametric inference based on BART under reasonable regularity conditions is still lacking.
Semiparametric inference based on BART is challenging because tree- and forest-based methods build on piecewise constant functions. While BART can capture additive structures and perform variable selection, the standard Donsker property, which requires the unknown function to have smoothness exceeding one-half of the covariate dimension, is difficult to justify. In the context of missing data and causal inference, we impose the standard identifying assumption that potential outcomes and the response indicator are independent, conditional on the observed covariates. Within this framework, the main objective of the paper is to develop a robust inference procedure for the mean response based on BART.
We propose a novel robust Bayesian procedure that is asymptotically valid without imposing the Donsker property or the no-bias condition. For estimation of the mean response in the missing data problem or the average treatment effect (ATE) in program evaluation, we explore the role of propensity score by making use of the double robust functional. By combining a preliminary estimator of the propensity score, the BART induced posterior for the conditional mean function, and the Bayesian bootstrap weights, we establish a semiparametric Bernstein-von Mises (BvM) Theorem, which implies inference on the mean response at a $\sqrt n$-rate and asymptotic efficiency in the semiparametric sense. When the Donsker property fails, however, an additional bias term appears in the posterior distribution in the Bernstein-von Mises Theorem. Interestingly, this bias term is of the exact same form as that obtained by BLY2022, who consider the prior correction approach of ray2020causal based on least favorable directions of semiparametric submodels.
To adjust for the bias in the BvM statement, we propose a posterior correction approach based on pilot estimators for the propensity score and the conditional mean functions. As a technical device, we rely on sample splitting for the pilot estimation, in line with the recent double machine learning (DML) literature ChernozhukovNeweySingh2020a. In doing so, our method differs from the proposal of yiu2023corrected, which places a prior on the propensity score or estimates it within a one-step posterior correction. In contrast to their BvM results (given in their Theorems 3 and 5), we do not impose Donsker conditions in our paper. This is possible because our additional posterior correction removes the bias term in non-Donsker regimes for the mean response. Such an adjustment is crucial for BART, which builds on non-smooth base learners and thereby only adapts to conditional mean functions that are Lipschitz smooth rockova2020posterior. As a result, the Donsker requirement is particularly stringent with multi-dimensional covariates. Our proposed robust BART (RoBART) achieves the BvM under more flexible smoothness conditions: the lack of smoothness of conditional mean functions can be compensated by high regularity of the propensity score.
Monte Carlo simulation results show that RoBART improves robustness over conventional BART while maintaining competitive credible interval lengths. Our method also shows finite-sample advantage over the one-step posterior correction. With increasing complexity of the underlying conditional mean and propensity score functions, one-step posterior correction is not sufficient to address the undercoverage of BART. By contrast, incorporating the additional debiasing term into the BART posterior allows our robust BART approach to achieve improved bias and coverage performance in complex model designs, without sacrificing estimation precision. We also showcase the practicality of our method by a real-data application. \\
Related literature. Our work is most closely related to active research areas. There has been considerable interest in developing semiparametric Bayesian inference that deviates from the plug-in principle. ray2020causal pioneered the use of data-dependent priors to weaken the regularity conditions on the propensity score, for which they coined the term single robustness. Maintaining the same prior construction, BLY2022 introduce an additional debiasing step that allows for double robustness. This is further extended to the popular difference-in-differences design by BLY2024. In this paper, we start from a different perspective by examining the one-step posterior of yiu2023corrected and establish its asymptotic equivalence with the prior adjusted method. More importantly, we conduct another posterior bias correction to remove the bias term $b_{0,\eta}$. These are designed to the relax the Donsker property of the conditional mean, which is required by yiu2023corrected, but it is not met for the BART with high-dimensional covariates and nonadditive functions.
As a general nonparametric function estimation, the BART achieves the adaptive near rate-minimax posterior contraction, which has been the major breakthrough, obtained by linero2018bart,rockova2020posterior,rockova2020bart,jeong2023art. Compared with the extensive results on the posterior convergence for BART, the corresponding inferential theory is much scarce. rockova2020bart established a semiparametric BvM theorem using BART type priors for a particular type of linear functional of the regression function, under the fixed design. Her result requires Donsker property of the conditional mean. In the presence of multi-dimensional covariates, she further requires the weight function to be uniform in this linear functional. For the Bayesian CART (for a single tree not the forest), castillo2021CART also establish general semiparametric BvM theorem within the Donsker regime.
The paper is structured as follows. Section (ref) describes the missing data model together with the identifying assumptions and presents the robust semiparametric Bayesian procedure. In Section (ref), we derive the BvM Theorem under high-level assumptions. Section (ref) illustrates these results for BART. Section (ref) provides numerical results via Monte Carlo simulations and an empirical illustration. Proofs of the main theoretical results are presented in Appendix (ref). Lemmas and auxiliary results are collected in Appendix (ref). Appendix (ref) extends the framework to inference on average treatment effects. Appendix (ref) provides additional implementation details of our methodology.\\
Notation. We adopt the standard empirical process notation as follows. For a function $h$ of a random vector $Z$ that follows distribution $P$, we let $P[h]=\int h(z)\,\mathrm{d} P(z),\mathbb{P}_n[h]=n^{-1}\sum_{i=1}^{n}h(Z_i)$, and $\mathbb{G}_n[h]=\sqrt n\left(\mathbb{P}_n-P\right)[h]$. For two sequences $\{a_n\}$ and $\{b_n\}$ of positive numbers, we write $a_n \lesssim b_n$ if $\limsup_{n\to\infty} (a_n / b_n)<\infty$, and $a_n \sim b_n$ if $a_n \lesssim b_n$ and $b_n \lesssim a_n$.
This section provides the main setup of our missing data framework and motivates the new Bayesian methodology. We provide the implementation of our robust Bayesian algorithm and discuss its connection and difference to the existing literature.
We consider a family of probability distributions $\{P_\eta:\eta\in\mathcal H\}$ for some parameter space $\mathcal H$. The (possibly infinite dimensional) parameter $\eta$ characterizes the probability model. An independent and identically distributed (i.i.d.) sample $\{(R_iY_i,R_i,X_i^\top)^\top\}_{i=1}^n$ is available, where the dependent variable $Y_i$ is observed only if $R_i=1$ and it is missing if $R_i=0$. We assume the outcomes $Y_i$ are missing at random (MAR); that is, the outcome $Y_i$ and the missingness indicator $R_i$ are conditionally independent given $X_i$.
Throughout the paper, the distribution of the observable vector $(RY,R,X^\top)^\top$ is completely modeled by three components: (i) the marginal distribution $\mu_{\eta}(x)$ of the $p$-dimensional covariates $X$; (ii) the propensity score $\pi_{\eta}(x)=P_{\eta}(R=1\mid X=x)$; and (iii) the conditional density function $f_{\eta}(\cdot \mid x)$ of the outcome $Y$ given $X=x$ and $R=1$. Consequently, the conditional mean function of the outcome given covariates is
Under the MAR assumption, the probability density function $p_{\eta}$ for $Z_i=(R_iY_i,R_i,X_i^\top)^\top$ can be written as
Let $\eta_0$ be the true value of the parameter and denote $P_0=P_{\eta_0}$, which corresponds to the frequentist distribution that generates the observed data. Accordingly, we denote $\mu_0$ as the true marginal covariates distribution, $\pi_0$ as the true propensity score, and $f_{0}(\cdot\mid x)$ as the true conditional density of $Y$ given $X=x$ and $R=1$.
The key parameter of interest is the mean of the outcome variable:
where $\mathbb{E}_0[\cdot]$ denotes the expectation under $P_0$. Under the MAR assumption, $\chi_0$ can be identified using the conditional mean function $m_0(x):=\mathbb{E}_0[R_iY_i\mid R_i=1,X_i=x]$. Specifically, we have
For our approach, it is useful to consider an alternative expression of $\chi_0$, augmented by a term depending on the inverse propensity score $\pi_0$ as follows
where the Riesz representer $\gamma_0$ is given by
which satisfies $\mathbb{E}_0[m_0(X)\gamma_0(R,X)]=\chi_0$; see chernozhukov2018double.
Let $W^{(n)}:=(W_{n1},\ldots, W_{nn})$ denote the Bayesian bootstrap weights by rubin1981bayesian, i.e., $W_{ni}=e_i/\sum_{j=1}^n e_j$ for $e_i \stackrel{iid}{\sim} \textup{Exp}(1)$, $i=1,\dots,n$. Referring to the expression (ref), the natural starting point for semiparametric Bayesian inference is
where $m_{\eta}(\cdot)$ is some stochastic function with a given prior and $\widehat{\gamma}(r,x)=r/\widehat{\pi}(x)$ is the estimated Riesz representer with a pilot estimator $\widehat\pi$ of $\pi_0$ computed over some auxiliary dataset. The posterior law of $\chi_{\eta}$ or its conditional distribution given the observed data $Z^{(n)}=(Z_1,\ldots, Z_n)$ is determined jointly by the posterior law of $m_{\eta}$ denoted by $\Pi_{m}(\cdot\mid Z^{(n)})$, as well as the distribution of random vector $W^{(n)}$ involving Bayesian bootstrap weights. For the latter part, its probability law is known and easy to simulate. In accordance with the study on the bootstrap law van1996empirical, we adopt the notation $\Pi_{W}(\cdot\mid Z^{(n)})$. While our general results are applicable to a broad class of nonparametric priors for $m_{\eta}$, we particularly analyzes prior modeling under more basic conditions based on the BART. The posterior draws of the vector $(m_\eta(X_1),\ldots,m_\eta(X_n))$ can be easily obtained using the R package $\mathsf{BART}$ sparapani2021nonparametric.
To motivate the approach from a Bayesian perspective, one can view the conditional distribution of ((ref)) as approximating the posterior of
where we use a degenerate posterior for the propensity score from the pilot estimation. When it comes to the prior for $P_{\eta}$, we consider the Dirichlet process prior with its base measure taken to be zero. It coincides with the Bayesian bootstrap process as in (ref). Note that the original notation in yiu2023corrected defines the one-step posterior of the mean functional $\chi_\eta$ given $Z^{(n)}$ as the conditional measure of
where the probability model $P_\eta$ is parameterized by $(m_\eta,\pi_\eta)$, and the augmented part $\tilde P$ denotes the Bayesian bootstrap law. The posterior law of $(m_\eta,\pi_\eta)$ is also assumed to be conditionally independent of the bootstrap law $\tilde{P}$. Both angles lead to the same approach because the Bayesian bootstrap law is known and does not depend on the unknown data generating process for the data. The main contribution of the current paper is that we show the one-step posterior still contains a bias term when the conditional mean function does not satisfy the Donsker property. This is particularly relevant for the BART with multi-dimensional covariates.
Our Bayesian approach connects to the least favorable direction as specified by the efficient influence function
for some Riesz representer $\gamma_\eta$ given by
A pilot estimator for the propensity score $\pi_0$ is denoted by $\widehat\pi$ based on an auxiliary sample, so that $\widehat{\gamma}(r,x)=r/\widehat \pi(x)$ is an estimator of the Riesz representer $\gamma_0$. We also make use of a pilot estimator $\widehat m$ for the conditional mean function $m_0$ in the debiasing step. The use of an auxiliary data for the estimation of unknown functional parameters simplifies the technical analysis and is common in the related Bayesian literature; see ray2020causal for propensity score adjusted priors in the case of missing data. It is also inspired by the recent literature on the DML chernozhukov2018double, which yield negligibly of certain smaller order terms. This technique also dates back to the early development in semiparametric estimation schick1986semi. In practice, we use the full data twice and do not split the sample, as we have not observed any over-fitting or loss of coverage thereby.
We are now in a position to describe our robust Bayesian procedure. The inputs to the algorithm are the observable data $\{(R_i Y_i, R_i, X_i^\top)^\top\}_{i=1}^n$, the pilot estimators $\widehat{\gamma}$ and $\widehat{m}$ as discussed above, and the number of posterior draws $S$. The posterior distribution of $(m_\eta(X_i))_{i=1}^n:=(m_\eta(X_1),\ldots, m_\eta(X_n))$ is generated by applying BART to estimate $m_{\eta}(\cdot)$ using the data $\{(R_iY_i, X_i^\top): R_i=1\}_{i=1}^n$.
Given the posterior simulation draws, the $100\cdot(1-\alpha)\%$ credible set $\mathcal{C}_n(\alpha)$ for the parameter of interest $\chi_0$ is computed by
where $q(a)$ denotes the $a$-quantile of $\{\check{\chi}_\eta^{s}:s=1,\ldots,S\}$. We take the Bayesian point estimator (the posterior mean) by averaging the simulation draws: $\overline{\chi}_{\eta}=S^{-1}\sum_{s=1}^S \check{\chi}_\eta^{s}$.
Given some pilot estimators $\widehat{m}$ for the conditional mean $m_0$ and $\widehat{\pi}$ for the propensity score $\pi_0$, the well-known doubly robust estimator takes the form:
This has been extensively studied in the frequentist literature, starting with the augmented inverse propensity score weighting (AIPW) by robins1994regression. Its extension to using pilot machine learning estimators combined with cross-fitting by chernozhukov2017double has generated considerable interest in the literature. The formulation in Step (ref) of our algorithm mimics the frequentist construction and it agrees with the one-step posterior proposed by yiu2023corrected.
Although our primary interest lies in establishing the BvM theorem for the BART, we first present the theory under high-level conditions. The robustness achieved through our additional debiasing step is thus applicable to other types of priors beyond the BART framework.
The one-step posterior of $\chi_{\eta}$ is determined jointly by the posterior of the conditional mean and the Bayesian bootstrap law. Recall that the conditional probability density function associated with the conditional mean function by $f_{\eta}(\cdot \mid x)$. For simplicity, we work with the conditional density class which only depends on the unknown conditional mean function. Such a class includes the standard Gaussian, Bernoulli, and Poisson outcomes, as discussed in BLY2022. Accordingly, we denote it by $f_{m_\eta}(\cdot \mid x)$ thereafter. The posterior distribution is formally given by
Regarding the centering point in the asymptotic normal approximation, we can consider any asymptotically efficient estimator $\widehat{\chi}$ with the following linear representation:
where $\widetilde{\chi}_0=\widetilde{\chi}_{\eta_0}$ is the efficient influence function given in (ref) under $\eta=\eta_0$. Denote its variance by $\textsc v_0=\mathbb{E}_0\left[ \widetilde{\chi}_0^2(Z)\right]$. We write $\mathcal{L}_{\Pi}(\sqrt{n}(\chi_\eta-\widehat{\chi})\mid Z^{(n)})$ for the marginal posterior distribution of $\sqrt{n}(\chi_\eta-\widehat{\chi})$.
Below, we consider some measurable sets $\mathcal H^m_n$ of functions $m_{\eta}$ such that $\Pi(m_{\eta}\in\mathcal{H}^m_n\mid Z)\to_{P_0} 1$. To abuse the notation for convenience, we also denote $\mathcal{H}_n=\{\eta:m_{\eta}\in\mathcal{H}_n^m\}$ where we index the conditional mean function $m_{\eta}$ by its subscript $\eta$. We introduce the notation $\|\phi\|_{ 2,\mu_0}:= \sqrt{\int \phi^2(x)\,\mathrm d\mu_0(x)}$ for all $\phi\in L^2(\mu_0):=\{\phi:\|\phi\|_{ 2,\mu_0}<\infty\}$, as well as the supremum norm $\|\cdot\|_\infty$. We introduce the notation for the conditional variance $\sigma_0^2(x)= \mathbb E_0[(RY-m_0(X))^2\mid R=1,X=x]$ and assume throughout the paper that $\sigma_0\in L^2(\mu_0)$. We denote the support of a random variable $V$ by $\mathcal V$.
Assumption (ref) imposes a missing-at-random (MAR) assumption and an overlap assumption, both of which are sufficient for identifying the conditional mean function $m_0(x)=\mathbb{E}_0[Y_i \mid X_i=x]$.
Assumption (ref) imposes sufficiently fast convergence rates for the estimators for the conditional mean function $m_0$ and the propensity score $\pi_0$. In practice, one can explore the recent proposals from ChernozhukovNeweySingh2020b and hirshberg2021augmented. The posterior concentration rate for the conditional mean can be derived by verifying high-level assumptions of ghosal2000rates.
Assumption (ref) (i) restricts the functional class $\mathcal{G}_n$ to form a $P_0$-Glivenko-Cantelli class; see Section 2.4 of van1996empirical. The moment conditions on the envelope functions follow from yiu2023corrected, which address the uniformity issues when showing the convergence of the conditional Laplace transform. Assumption (ref) (ii) imposes a new stochastic equicontinuity condition on the product structure involving $\widehat\gamma$ and $m_\eta$, which sets us apart from the existing literature. In contrast to this assumption, ray2020causal and yiu2023corrected both require a stochastic equicontinuity condition $\sup_{\eta\in\mathcal{H}^m_n}\left|\mathbb{G}_n\left[m_\eta-m_0\right]\right|=o_{P_0}(1)$\footnote{Because a nonparametric prior is also assigned to the propensity score, yiu2023corrected further requires the propensity score function to belong to the Donsker class.}. As demonstrated by BLY2022, this is weaker than directly imposing the Donsker property on $(m_{\eta}-m_0)$. For the H\"older smooth class, BLY2022 show that the low-order differentiability of the conditional mean can be compensated by exploring the high-order smoothness of the propensity score, and vice versa. Note that $\sup_{\eta\in\mathcal{H}^m_n}\left|\mathbb{G}_n\left[m_\eta-m_0\right]\right|$ diverges for the non-Donsker class, by the Sudakov's inequality van1996empirical.\\
We now present a semiparametric Bernstein–von Mises theorem, which establishes asymptotic normality of the posterior distribution, modulo a bias term. This asymptotic equivalence result is established using the so called bounded Lipschitz distance on probability distributions on $\mathbb R$, see (see Chapter 11 of dudley2002). As the one-step posterior is built from mimicking the frequentist double robust estimand, it is surprising at first sight that this does not fully remove the bias term, hence the technical proof of the BvM theorems in yiu2023corrected still relies on the Donsker property. Note that one crucial difference to the frequentist method is that we need to characterize the weak convergence for the entire posterior distribution (not just its center as the point estimator), there is this crucial bias term that we identify. This shows the distinctive feature of deriving the frequentist validity of semiparametric Bayesian inference.
The previous result shows that, under double robust smoothness conditions, the posterior of the one-step corrected mean only satisfies the BvM result within a bias term $b_{0,\eta}$. We emphasize that the Bayesian procedure that achieves the BvM equivalence in Theorem (ref) is not feasible, because it depends on the term $b_{0,\eta}$, which is a function of the unknown conditional mean $m_0$. This bias term vanishes, if we impose additional smoothness restriction on the conditional mean function $m$ satisfying Donsker property.
Aimed at relaxing the Donsker property, our objective is to maintain double robust smoothness conditions while considering pilot estimators for the unknown functional parameters in $b_{0,\eta}$. We explicitly correct the posterior distribution, following our proposed methodology in Algorithm (ref). Recall the definition of the bias correction term $\widehat{b}_{\eta}$ given in Algorithm (ref):
The next theorem is an immediate consequence of BLY2022.
As importance consequences of the BvM theorems, we hence provide the frequentist validity of the Bayesian credible set built from the corrected posterior. For any $\alpha\in(0,1)$ and Bayesian credible set $\mathcal{C}_n(\alpha)$ for $\chi_\eta-\widehat b_{0,\eta}$, then we have $ P_0\big(\chi_0\in \mathcal{C}_n(\alpha)\big) \to 1-\alpha.$
We turn to our main task of establishing a new BvM theorem where the conditional mean function is modeled by the BART. The BART prior expresses the unknown function as a sum of many binary regression trees. Each individual tree consists of a set of internal decision nodes which define a partition of the covariate space, as well as a set of terminal nodes or leaves. Those partitions are generated by recursively applying some binary split rules of the covariate space as the form $\{x_j \leq \tau \}$ versus $\{x_j>\tau \}$ where $x_j$ signifies that $j$-th coordinated is chosen for the splitting and $\tau$ is the split point. The good approximation property of tree-based learners relies on finding a fine partition scheme that divides the data into more homogeneous groups and learning a piecewise constant function on the partition. For this purpose, we introduce the notion about the split-net from rockova2020posterior. Given an increasing integer sequence $b_n$, a split-net $\mathcal{X}=\{x_j\in[0,1]^p,j=1,\ldots,b_n \}$ is a discrete collection of $b_n$ points $x_j$ at which possible splits occur along coordinates. For a set $\Omega\subset\mathbb{R}^p$, we denote the $j$--th projection mapping of $\Omega$ by $[\Omega]_j=\{x_j\in \mathbb{R}:(x_1,\ldots,x_p)^\top\in \Omega \}$.
The construction of each individual tree starts with some root node in $[0,1]^p$ at depth $l=0$, where the depth of a node means the number of nodes along the path from the root node down to that node. Any binary tree can be characterized by the (i) the probability that a node at depth $l$ is nonterminal, (ii) the distribution on the splitting variable assignments at each splitting node, and (iii) the rules with which the splits are made. Referring to point (i), each node at depth $l\in\{0,1,2,\ldots\}$ is split (hence nonterminal) with some prior probability depending on the depth. When it comes to point (ii), we follow linero2018dart,linero2018bart using a sparse Dirichlet prior. This approach chooses a splitting covariate $j$ from a proportion vector $\vartheta=(\vartheta_1,\ldots,\vartheta_p)^\top$ belonging to the $p$-dimensional simplex. Finally for point (iii), if a node corresponding to a box $\Omega$ is split, a splitting coordinate is drawn from the proportion vector and a split-point $\tau_j$ is chosen randomly from $[\mathcal{X}]_j\bigcap\textsf{int}([\Omega]_j)$ for a given split-net $\mathcal{X}$. The procedure ends until all nodes become terminal. We further denote the cardinality of the split points in the split-net $\mathcal{X}$ projected onto the $j$-th coordinate by $b_j(\mathcal{X})$ for $j=1,\ldots,p$.
Throughout the paper, we pick a fixed number of trees $T$ and assume an independent product prior for the tree ensemble, $\Pi(\mathcal{E})=\prod_{t=1}^T \Pi(\mathcal{T}^t)$. chipman2010bart recommend taking $T=200$ based on extensive numerical evidence. For a given split-net $\mathcal{X}$ and for each $1\leq t\leq T$, we denote with $\mathcal{T}^t=(\Omega^t_1,\ldots,\Omega^t_{K^t})$ a $\mathcal{X}$-tree partition of size $K^t$ and with step-heights as $\bm{\beta}^t=(\beta_1^t,\ldots,\beta_{K^t}^t)\in\mathbb{R}^{K^t}$. An additive tree-based learner is fully described by a tree ensemble $\mathcal{E}=\{\mathcal{T}^1,\ldots,\mathcal{T}^T \}$ and terminal node parameters $\mathcal{B}=(\bm{\beta}^{1},\ldots,\bm{\beta}^{T})^\top\in\mathbb{R}^{\sum_{t=1}^TK^t}$ as follows
Given an ensemble $\mathcal{E}$ of trees, we denote $\mathcal{M}_{\mathcal{E}}:=\{ m_{\mathcal{E},\mathcal{B}}(x):\mathcal{B}\in\mathbb{R}^{\sum_{t=1}^TK^t} \}$. If $\mathcal{E}$ consists of a single tree $\mathcal{T}$, we denote $\mathcal{M}_{\mathcal{E}}$ by $\mathcal{M}_{\mathcal{T}}$. Binary decision trees or forests offer good approximation properties to the class of H\"older continuous class, which is indexed by an exponent $\alpha\in(0,1]$. This parameter $\alpha$ does not exceed one, which is standard in the literature studying the piecewise constant estimators and priors. Another notable feature of BART is that one can explore the sparsity of relevant regressors and perform variable selection. The true conditional mean function is assumed to belong to the following class of functions. Below we introduce the norm $\Vert \phi\Vert_{\mathcal{H}^{\alpha}}:= \sup_{x,y\in[0,1]^p}|\phi(x)-\phi(y)|/\Vert x-y\Vert_2^{\alpha}$, for functions $\phi:[0,1]^p\mapsto \mathbb{R}$, with some $\alpha>0$, where $\|\cdot\|_2$ denotes the Euclidean norm.
For the set $\mathcal{S}_0$ with its cardinality $|\mathcal{S}_0|=q_0<p$, we assume $m_0\in \mathcal{F}_p(\alpha,\mathcal{S}_0)$. In addition, we also explore the case where $m_0$ is additively separable in terms of low-dimensional covariates. Consider the following additive class with $T_0$ components
where $\bm{\alpha}:=(\alpha^t)_{t=1}^{T_0}$ and $\bm{\mathcal{S}}:=(\mathcal{S}^t)_{t=1}^{T_0}$. We denote $q_0^t=|\mathcal{S}^t_0|$ for the subset $\mathcal{S}^t_0\subset\{1,\ldots,p\}$. We define
related to the rate of posterior contraction for $m_0\in \mathcal{F}_p(\alpha,\mathcal{S})$. When considering the additive case $m_0\in \mathcal{F}^{add}_p(\bm{\alpha},\bm{\mathcal{S}})$, we define $\varepsilon^{add}_{n}=\sqrt{\sum_{t=1}^{T_0}\varepsilon_{n,t}^2}$, where $\varepsilon_{n,t}:=n^{-\alpha^t/(2\alpha^t+q^t_0)}\sqrt{\log n}$.
The case where $m_0\in \mathcal{F}_p(\alpha,\mathcal{S}_0)$ provides a clean illustration of our theory, which shows that our additional debiasing step is essential if $q_0>1$. The class $\mathcal{F}_p(\alpha,\mathcal{S}_0)$ fails to satisfy the Donsker property whenever $\alpha<q_0/2$. A subtle consequence, is that the size of the resulting sieve set used to approximate the H\"older continuous function becomes too large, so that the standard maximal inequality via the entropy bound fail to deliver the stochastic equicontinuity. The same issue occurs to the additive class whenever $\alpha^t<q^t_0/2$ for any $1\leq t\leq T_0$.
Next, we list assumptions of covariates, error term and the prior specifications, followed by some remark.
Theorem (ref) establishes a BvM result for BART without requiring the Donsker property for either the propensity score or the conditional mean function. This result follows from verifying the high-level assumption from the previous section. The proof is based on constructing a sieve set that receives posterior mass with probability approaching one. This sieve set is defined by piecewise constant functions over increasingly fine grids and with a growing number of covariates. Due to the discontinuity of such piecewise constant functions, the complexity with increasing number of covariates exceeds the threshold of the Donsker class. Hence, the standard stochastic equicontinuity as imposed in the semiparametric Bayesian literature is not satisfied.
In this section, we first examine the finite-sample performance of BART inference for the mean response $\chi_0=\mathbb{E}_0[Y_i]$ in missing data models. To illustrate the practical relevance of the theoretical discussion in Remark (ref), we consider simulation designs that include interactions. As Appendix (ref) extends the proposed robust BART inference to the average treatment effect (ATE) framework, we then apply it to the well-known National Health and Nutrition Examination Survey data to revisit the average effect of participation in meal programs on students’ body mass index.
The data-generating process for i.i.d. observations is specified as follows. For each unit $i=1,\ldots, n$, we generate the observed variables $(R_iY_{i}, R_i, X_{i}^\top)$ by $X_i = (X_{i1},\dots, X_{i5})^\top$ where $X_{i1}, X_{i2}, X_{i3}\sim N(0,1)$, $X_{i4}\sim \text{Bernoulli}(0.5)$, $X_{i5}$ is a categorical variable taking values $\{1, 2, 3\}$ with equal probability. The distributions of $R_i$ and $Y_i$ are given by
where $\Psi(t)=1/(1+e^{-t})$. We analyze the finite-sample effects of varying the complexity of the conditional mean function $m$ and the propensity scores for sample sizes $n \in \{125, 250, 500\}$. Throughout our simulations, the number of Monte Carlo replication is set to 1000.
We consider four designs for the function $e(\cdot)$ that determines the propensity score by $\pi(x)=\Psi(e(x))$ and the conditional mean function $m(\cdot)$. In Design I, $e(x)$ includes an interaction term $x_1x_3$, whereas $m(\cdot)$ contains no interaction terms. In Designs II–IV, both $e(\cdot)$ and $m(\cdot)$ include interaction terms, with increasing complexity (i.e., more interaction terms) across designs:
where $h(x_5)= 2\mathbbm 1\left\{x_5=1\right\}-\mathbbm 1\left\{x_5=2\right\}-\mathbbm 1\left\{x_5=3\right\}/2$ with $\mathbbm1\{\cdot\}$ denoting the indicator function.
We evaluate three inference methods. Standard BART obtains the posterior of the conditional mean function $m_{\eta}(\cdot)$ using BART, and then averages over the covariates using Bayesian bootstrap weights. Implementation relies on the R package $\mathsf{BART}$ sparapani2021nonparametric: the posterior of $m_{\eta}(\cdot)$ is computed using the function $\mathsf{gbart}$ with argument $\mathsf{type= wbart}$. We draw $2000$ posterior samples after discarding a burn-in of $500$, with the number of trees set to $T=200$. Default priors from the package are used, see Appendix (ref) for details. One-step BART applies the posterior correction of yiu2023corrected, which uses posteriors of both $m(\cdot)$ and $\pi(\cdot)$ obtained via BART or its logistic variant.\footnote{We implement logistic BART by calling $\mathsf{gbart}$ with $\mathsf{type= lbart}$.} RoBART, the robust BART method described in Algorithm 1, combines the BART posterior of $m_{\eta}(\cdot)$, a plug-in estimator for $\pi(\cdot)$, and the debiasing term $\widehat{b}_{\eta}$ in ((ref)). We consider two estimators for propensity score: Logit, a quadratic logistic regression; and SL, a SuperLearner combining logistic regression and a generalized additive model.\footnote{We use the R package $\mathsf{SuperLearner}$ with the library of prediction algorithms including $\mathsf{SL.glm}$ and $\mathsf{SL.gam}$ for implementation.} Both include pairwise interactions of covariates.
Table (ref) presents the bias of the posterior mean, coverage probability (CP) and the average length (CIL) of the $95\%$ credible interval formed by the (corrected) posteriors. In Design I, where $m(\cdot)$ is linear, all methods perform well including standard BART. In Design II, where interaction terms appear in both $e(\cdot)$ and $m(\cdot)$, BART tends to undercover, while one-step BART and RoBART restore CP close to nominal levels.\footnote{This complements the simulations in yiu2023corrected, which show improved coverage for one-step BART in a design where both $\pi(\cdot)$ and $m(\cdot)$ are one-dimensional but discontinuous, see their Table 2.} Between the two, RoBART exhibits smaller bias, especially in small samples. In Designs III and IV, as more interaction terms appear in $e(\cdot)$ and/or $m(\cdot)$, one-step BART produces larger bias and undercoverage, whereas RoBART continues to yield small bias and improved coverage. These results align with our theory. As the additive components of $\pi(\cdot)$ or $m(\cdot)$ depend on more than one covariate (see Remark (ref)), the Donsker property is not guaranteed. Evidently, our RoBART demonstrates robust performance in such complex designs. In addition to its improved bias and coverage performance, RoBART (with a logit estimator for the propensity score) yields shorter credible intervals than one-step BART for all designs and sample sizes, while RoBART (with a SuperLearner for the propensity score) does so in 9 out of 12 cases.
We apply the BART inference methods to a subsample of data from the National Health and Nutrition Examination Survey (NHANES) 2007–2008, previously analyzed by chan2016globally. The parameter of interest is the average treatment effect (ATE), defined as $\mathbb{E}[Y_i(1)-Y_i(0)]$, of participating in school meal programs ($D$) on body mass index (BMI) for children and youths aged 4--17 years ($Y$). The covariates $X$ include eleven variables: child age, child gender, race dummies (Black and Hispanic), an indicator for family income above $200\%$ of the federal poverty level, indicators for participation in the special supplemental nutrition program and in the food stamp program, an indicator for childhood food security, an insurance coverage dummy, and the age and gender of the survey respondent (an adult in the family). The sample size is $2330$.
We adapt the BART inference methods from the simulation study to the ATE setting (see Appendix (ref) for the robust BART (RoBART) method for ATE), and present the results in Table (ref). Table (ref) shows that all ATE estimates of school meal programs on BMI are small relatively to the sample average BMI $20.11$ (with the standard deviation $5.42$). All $95\%$ credible intervals include zero, suggesting that the average effect may go either direction when uncertainty is taken into account. One-step BART produces more dispersed posterior and thus more volatile point estimates than other methods. The estimates tend to increase when observations with propensity score near the boundary are discarded. We compare the BART-based estimates in Table (ref) with the frequentist results reported in Table 1 of chan2016globally using the full sample ($\bar n=2330$): Horvitz-Thompson (HT) estimate of $-1.48$ with a $95\%$ confidence interval $[-2.50,-0.46]$, the inverse probability weighting (IPW) estimate $-0.41$ with confidence interval $[-0.62, 0.34]$, and their preferred calibration estimate of $-0.04$ with the confidence interval $[-0.48, 0.40]$.\footnote{We cite the calibration estimator with exponential tilting as representative; other calibration estimates in chan2016globally are similar in magnitude and interval length.} While point estimates differ in sign across methods, all but the HT estimator indicate that participation in school meal programs had no statistically significant effect on children’s and youths’ BMI. In particular, both our proposed RoBART and the calibration estimator of chan2016globally yield small point estimates when using the full sample. And the credible intervals generated by RoBART are similar in magnitude to the confidence intervals of the calibration estimator.