Bayesian Semi-supervised Inference via a Debiased Modeling Approach

Bayesian Semi-supervised Inference via a Debiased Modeling Approach

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

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

{\bf Keywords:} Robustness and efficiency; Bayesian semi-parametric inference; Debiasing; Sample splitting and cross-fitting; Bernstein--von Mises theorem.

Introduction and overview of contributions

Semi-supervised (SS) learning has emerged as an exciting and active research area in statistics and machine learning in recent years. A typical SS setting involves two types of data sets: (i) a small or moderate sized labeled (or supervised) data $\mathcal{L}$ {with} observations for both an outcome (or label) $Y$ and a set of predictors $\mathbf{X}$, and (ii) a much larger sized unlabeled (or unsupervised) data $\mathcal{U}$ containing observations only for $\mathbf{X}$. SS settings arise naturally when the outcome is difficult or costly to obtain, but observations for the predictors are plenty and easy to access. Typically, this scenario occurs in many modern big-data problems involving large (electronic) databases, such as speech recognition, text mining, and more recently, biomedical applications like electronic health records books/mit/06/CSZ2006, zhu05survey, kohane2011using, chakrabortty2018efficient. In a standard SS setup, one of the primary statistical goals is to investigate whether and how parameter estimation and accuracy of inference can be improved by making use of the unlabeled data $\mathcal{U}$, unlike supervised methods, which use only the labeled data $\mathcal{L}$ and completely ignore $\mathcal{U}$. SS inference in this spirit has been studied in the recent frequentist literature for various problems, including mean estimation zhang2019semi,zhang2022high and linear regression chakrabortty2018efficient,azriel2022semi, among others. However, Bayesian approaches for SS inference are largely lacking in the literature to the best of our knowledge.

We propose a Bayesian debiased modeling and inference (BDMI) procedure for estimating the {\it population mean} $\theta_0:= \mathbb{E}(Y)$ of $Y$ under the SS setting, as a prototypical example. {A fundamental idea behind BDMI is to} carefully model certain {\it summary statistics} of the data in a {\it targeted} manner, rather than specifying a probability model for the raw data itself, along with developing and exploiting a novel {\it Bayesian notion of debiasing of nuisance parameters} (that are inherently involved in the procedure). Most existing SS approaches for estimating $\theta_0$ (or similar parameters/functionals of the distribution of $Y$) naturally require estimation of the possibly high dimensional {\it regression function} $m_0(\mathbf{X}):= \mathbb{E}(Y|\mathbf{X})$ to exploit $\mathcal{U}$ chakrabortty2018efficient, zhang2019semi, TonyCai2018SemisupervisedIF, zhang2022high. $m_0(\cdot)$ therefore acts as a {\it nuisance function} here, that is needed (for exploiting $\mathcal{U}$) but is not of primary interest. In general, the presence of such a nuisance parameter and its own estimation bias can drastically affect the final estimator's asymptotic behavior in the first order. In recent years, a popular frequentist debiasing procedure called double machine learning (DML) based on Neyman orthogonalization has been developed to rectify the impact of bias in learning a nuisance parameter chernozhukov2018double. A key contribution of this work is to develop a {\it Bayesian analogue} of such debiasing procedures, that ensures {\it robust, efficient} and {\it nuisance-insensitive} Bayesian inference for $\theta_0$ (the {\it target}) while {\it allowing for slow/inefficient (or even inconsistent) learning of} $m_0$.

BDMI encapsulates a new principle of {\it disentangling} the nuisance parameter that is amenable to {\it Bayesian} modeling and inference. It crucially relies on a debiased representation (Section (ref)) of $\theta_0$ in terms of $m_0$ (specifically, its estimator or a posterior sample) that simultaneously exploits $\mathcal{U}$ and also captures the nuisance bias incurred. Exploiting this representation, we then propose to model carefully chosen summary statistics of the data (see Section (ref)). Modeling summary statistics of the data has been sporadically considered in the Bayesian literature for estimation and hypothesis testing pratt1965bayesian, savage1969nonparametric, doksum1990consistent, clarke1995posterior,johnson2005bayes,lewis2021 as well as in likelihood-free inference methods like Approximate Bayesian Computation (ABC) marjoram2003markov, fearnhead2012constructing, drovandi2015. In the present setting, the summary statistics are exploited to: (i) carefully pinpoint the target and the bias induced from the nuisance, and (ii) learn them jointly by constructing a robust working likelihood (that can be justified under mild assumptions on the data generating mechanism) which can then be combined with default prior distributions on the model parameters to arrive at a posterior distribution. Further, a key feature of our approach is the careful usage of {\it sample-splitting} and {\it cross-fitting (CF)} chernozhukov2018double, newey2018cross -- {\it not} just as a technical artifact (as is common in the frequentist literature) but as an {\it integral component} of the debiasing process itself. It helps create {\it independent} sub-folds of the entire data that crucially enable the disentangling of the nuisance estimation process from the summary statistics modeling process. Further, to ensure usage of the full data overall, we use CF by rotating the roles of the splits and using each sub-fold in turn, and thereafter {\it aggregating} the posteriors from all sub-folds using a consensus Monte Carlo type approach scott2022bayes. It is worth mentioning that, while commonplace in the modern frequentist literature on semi-parametric inference, handling sample splitting (and CF) under a Bayesian framework is more challenging since it requires combining {\it distributions} (posteriors) and not just point estimators. Our final CF-based version of BDMI is given in Section (ref) and summarized in Algorithm (ref).

We show through our theoretical results in Section (ref) that the marginal posterior distribution $\Pi_{\boldsymbol \theta}$ for $\theta$ from BDMI inherits a {\it Bernstein--von Mises (BvM)-type limiting behavior} van2000asymptotic with an asymptotically Gaussian shape, and contracts {\it always} around the true $\theta_0$ at a parametric $n^{-1/2}$ rate ($n$ being the size of $\mathcal{L}$) and with a spread {\it tighter} than the supervised counterpart -- all holding {\it irrespective} of the choice/method used to obtain the nuisance posterior ($\Pi_{\mathbf{m}}$) for learning $m_0$. Further, $\Pi_{\boldsymbol \theta}$'s first order variability is {\it unaffected} by that of $\Pi_{\mathbf{m}}$ and is of the correct $n^{-1/2}$ rate {\it even if} the contraction rate of $\Pi_{\mathbf{m}}$ is arbitrarily {\it slow} or if it is even {\it misspecified} (i.e., does not contract around the true $m_0$). This makes BDMI first-order insensitive chernozhukov2018double to the nuisance estimation. Most importantly, from an SS inference perspective, $\Pi_{\boldsymbol \theta}$ (and its posterior mean) provably possess the desirable properties of global robustness and efficiency improvement: we show (i) the symmetric Bayesian credible intervals (CIs) from $\Pi_{\boldsymbol \theta}$ possess asymptotically correct frequentist coverage and sizes (of order $n^{-1/2}$) guaranteed to be {tighter} than their supervised counterpart; and (ii) the posterior mean is always $\sqrt{n}$-consistent, asymptotically Normal and more efficient or at least as efficient as the supervised estimator. Furthermore, when $\Pi_{\mathbf{m}}$ is correctly specified (with arbitrary contraction rate), $\Pi_{\boldsymbol \theta}$ and its posterior mean attain {\it optimal} efficiency, with variance matching the {\it semi-parametric efficiency bound}. All our claims above are validated through extensive simulations as well as a real data application in Section (ref). It is also worth noting that BDMI is {\it computationally scalable,} with all ingredient posteriors (from each fold) in $\Pi_{\boldsymbol \theta}$ being convolutions of $t$-distributions (hence {easy to sample} from). To our knowledge, BDMI is the {\it first work} on Bayesian inference (with provable guarantees) in SS settings.

Aside from SS inference itself, this work also contributes more generally to the growing literature on {\it Bayesian semi-parametric inference} in modern big-data settings. The SS setting has a distinct semi-parametric flavor, with $m_0(\cdot)$ being the (potentially high dimensional) nuisance parameter and functionals like $\theta_0$ being the target. There is a growing literature on frequentist properties of Bayesian semi-parametric inference procedures; see, e.g., bickel2012semiparametric, rivoirard2012bernstein,castillo2015bernstein, norets2015bayesian, ray2019debiased; where the quantity of interest is the marginal posterior of the parameter of interest obtained upon marginalizing out the nuisance parameter. Under delicate conditions on the prior distribution of the nuisance parameter, BvM results have been established for the parameter of interest in some of these works. Moreover, there have been some recent developments in the Bayesian semi-parametric literature (primarily for missing data or causal inference problems) aimed at alleviating bias arising from the nuisance estimation with slow rates ray2020semiparametric,luo2023semiparametric,breunig2022double,yiu2023semiparametric. Most of these are based on careful prior selection/modification, or tailored posterior updating, to mimic the flavors of their frequentist counterparts. BDMI adds to this literature by considering a different perspective and a principled approach to mitigate the bias of nuisance parameters. Another key feature of the approach is that it leaves the nuisance estimation method {\it entirely} to the user, and the nuisance posterior (or prior) does {\it not} require any form of adjustment or updating. While proposed albeit under the auspices of the SS inference problem, we believe the fundamental ideas of BDMI -- Bayesian debiasing and targeted modeling via summary statistics -- will also apply more generally to other Bayesian semi-parametric inference problems.

The rest of the article is organized as follows. We discuss the problem setup and some key preliminaries in Section (ref). Our proposed methodology is presented in Section (ref), with its various facets distributed across Sections (ref)--(ref). The theoretical properties of our method, including our main results (Theorems (ref)--(ref)), are presented in Section (ref), along with an alternative hierarchical version of our method and its theoretical properties discussed in Section (ref). Finally, extensive simulation studies and real data analysis are presented in Section (ref) to illustrate its empirical performance, followed by a concluding discussion in Section (ref). All technical materials, including the proofs of all the main theoretical results, along with supporting lemmas and their proofs, as well as additional numerical results and methodological discussions that could not be accommodated in the main paper, are collected in the \hyperref[sec:supplementary]{Supplementary Material} (Sections (ref)--(ref)).

The problem setup and key preliminary ideas

Let $Y \in \mathbb{R}$ be the outcome variable, $\mathbf{X} \in \mathbb{R}^p$ be the covariate (or predictor) vector, and $\mathbb{P}_{\mathbf{Z}} \equiv \mathbb{P}_{Y \mid \mathbf{X}} \otimes \mathbb{P}_{\mathbf{X}}$ be the unknown joint distribution of $\mathbf{Z} := (Y, \mathbf{X}')'$, where $\mathbb{P}_{Y \mid \mathbf{X}}$ and $\mathbb{P}_{\mathbf{X}}$ denote the conditional distribution of $Y\mid \mathbf{X}$ and the marginal distribution of $\mathbf{X}$, respectively. The {\it available data} under the SS setting is denoted as: $\mathcal{D} := \mathcal{L} \hspace{0.7mm} \cup \hspace{0.7mm} \mathcal{U}$, with $\mathcal{L} := \{ \mathbf{Z}_i \equiv (Y_i, \mathbf{X}_i')' : i = 1, \dots , n \}$ being the labeled data containing $n$ independent and identically distributed (i.i.d.) samples of $\mathbf{Z} \sim \mathbb{P}_{\mathbf{Z}}$, and $\mathcal{U} := \{\mathbf{X}_i : i= n+1, \dots, n+N \}$ being the unlabeled data containing $N$ i.i.d. samples of $\mathbf{X} \sim \mathbb{P}_{\mathbf{X}}$, and $\mathcal{L}$ and $\mathcal{U}$ are independent, denoted as $\mathcal{L} \perp \!\!\! \perp \mathcal{U}$.

\remark Assumption (ref) is fairly standard in the SS inference literature books/mit/06/CSZ2006, kawakita2013semisupervised. The condition (i) encodes a key (and {\it unique}) feature of SS settings, allowing for disproportionate sizes of $\mathcal{L}$ and $\mathcal{U}$. For example, while the size of $\mathcal{L}$ may be of the order of hundreds, the size of $\mathcal{U}$ could be of the order of tens of thousands. Further, since the outcome $Y$ is missing in $\mathcal{U}$, one can view SS inference as a missing data problem by assuming $Y$ is `missing completely at random' tsiatis2006semiparametric. However, since $\lim_{n, N \to \infty} n/N \to c = 0$ is allowed, it naturally violates the positivity assumption (on the proportion of $Y$ observed) standard in the missing data literature tsiatis2006semiparametric, and makes the SS setting fundamentally {\it different} and more challenging (due to {\it non-standard asymptotics}) from the missing data setup. The condition (ii) asserts that the underlying distributions of $\mathcal{L}$ and $\mathcal{U}$ are the same, which is standard and often implicit in the SS inference literature kawakita2013semisupervised, chakrabortty2018efficient, zhang2019semi, zhang2022high, along with a mild moment assumption on $\mathbf{Z}$ to ensure $\mathbb{E}(Y\mid\mathbf{X})$ and $\mathbb{E}(Y)$ exist. Finally, we clarify that we allow {\it high dimensional} settings throughout ($p$ can diverge with $n$).

Preliminaries: Notational conventions and the supervised approach

We use the following {\it notational conventions} throughout the paper. Let $\mathbb{E}(\cdot) \equiv \mathbb{E}_{\mathbf{Z}}(\cdot)$, $\mathbb{E}_{Y\mid\mathbf{X}}(\cdot)$ and $\mathbb{E}_{\mathbf{X}}(\cdot)$ denote expectations under the distributions $\mathbb{P} \equiv \mathbb{P}_{\mathbf{Z}}$, $\mathbb{P}_{Y\mid\mathbf{X}}$ and $\mathbb{P}_{\mathbf{X}}$, respectively. For any dataset/collection (or its subset/functions) $\cal{C}$ on $\mathbf{Z}$, let $\mathbb{E}_{\cal{C}}(\cdot)$ and $\mathbb{P}_{\cal{C}}(\cdot)$ denote expectations and probability under the joint distribution of $\cal{C}$. Let $W$ be a generic random variable (or vector) with an underlying probability distribution $\mathbb{P}_{W}$, and let $f$ be any measurable $\mathbb{R}$-valued deterministic function of $W$. Then, the expectation of $f(W)$ is defined as: $\mathbb{E}_{W}\{f(W)\} \equiv \mathbb{E}_{W \sim \mathbb{P}_W}\{f(W)\} := \int f(w) \mathrm{d} \mathbb{P}_{W}(w)$, whenever the Lebesgue integral exists. Further, for any \( d \geq 1 \), let \(\mathbb{L}_d(\mathbb{P}_{\mathbf{Z}})\) and \(\mathbb{L}_d(\mathbb{P}_{\mathbf{X}})\) denote the spaces of all \(\mathbb{R}\)-valued measurable functions \(g\) of $\mathbf{Z}$, and \(h\) of $\mathbf{X}$, such that \( \|g(\mathbf{Z})\|^d_{\mathbb{L}_d( \mathbb{P}_{\mathbf{Z}})} := \mathbb{E}_{\mathbf{Z}}\{| g(\mathbf{Z}) |^d \} < \infty\) and \( \|h(\mathbf{X})\|^d_{\mathbb{L}_d( \mathbb{P}_{\mathbf{X}})} := \mathbb{E}_{\mathbf{X}}\{| h(\mathbf{X}) |^d \} < \infty\), respectively. Let $\mathcal{N}(\mu, \sigma^2)$ denote the Normal (Gaussian) distribution with mean $\mu$ and variance $\sigma^2$, and $t_\nu(\mu, c^2)$ denote the $t$-distribution with degrees of freedom $\nu > 0$, center $\mu$ and scale $c$. We also use $\mathcal{N}(x; \mu, \sigma^2)$ and $t_\nu(x; \mu, c^2)$ to denote their respective probability density functions (pdfs) evaluated at $x \in \mathbb{R}$. For given probability measures $P$ and $Q$ on a measurable space $(\Omega, \mathcal{F})$, the total variation (TV) distance between $P$ and $Q$ is $\| P - Q \|_{\mathrm{TV}}:= \sup_{B \in \mathcal{F}} |P(B) - Q(B)|$. For a sequence $b_n > 0$ and a sequence of random variables $X_n$, we say $X_n = o_{\mathbb{P}}(b_n)$ if and only if (iff) $|X_n|/b_n \overset{\mathbb{P}}\to 0$ as $n \to \infty$. If $X_n \overset{\mathbb{P}}\to 0$, we write $X_n = o_{\mathbb{P}}(1)$. Similarly, a sequence of random variables $W_n = O_{\mathbb{P}}(b_n)$ iff for any $\varepsilon > 0$, there exist $B_\varepsilon > 0$ and $n_\varepsilon$ such that $\mathbb{P}(|W_n| \leq B_{\varepsilon} \hspace{0.5mm} b_n) > 1 - \varepsilon$ for all $n \geq n_{\varepsilon}$. Furthermore, $W_n = o_{\mathbb{P}}(1)$ iff for some sequence $b_n \to 0$, $W_n = O_{\mathbb{P}}(b_n)$. Lastly, for $\psi_0 \equiv \psi_0(\mathbb{P})$ denoting any functional of interest for any distribution $\mathbb{P}$, we let $\psi$ represent the corresponding {\it random variable} (or vector, function, etc., as applicable) in a Bayesian framework, and denote its posterior distribution by $\Pi_{\boldsymbol{\psi}}$. {\it This convention is used consistently, without mention, throughout the paper.}

Before discussing any SS approaches, we first introduce the standard {\it supervised} Bayesian approach for estimating $\theta_0$ using $\mathcal{L}$ only, to set a benchmark. In the supervised setting, one can adopt a Bayesian framework by modeling $\mathcal{L}$ (i.e., the $Y_i$'s $\in \mathcal{L}$) with a working Gaussian likelihood with mean $\theta$ and variance $\sigma^2$, combined with a joint prior on $(\theta, \sigma^2)$. This yields a marginal {\it posterior} $\Pi_{\sup}$ for $\theta$ which, under mild regularity conditions on the prior, satisfies a BvM result van2000asymptotic: $\Pi_{\sup} \approx \mathcal{N}(\widehat \theta_{\sup}, \sigma^2_{Y}/n)$ as $n \to \infty$, where $\widehat \theta_{\sup}:= \overline Y \equiv n^{-1}\sum_{i = 1}^n Y_i$ and $\sigma^2_Y:= \mathrm{Var}(Y)$. Thus, $\Pi_{\sup}$ yields $\widehat \theta_{\sup} \equiv \overline Y$ as a natural (supervised) {\it point estimator} of $\theta_0$, as well as CIs of sizes $\propto \sigma_Y/\sqrt{n}$. Further, $\sigma^2_Y$ is the {\it best} achievable variance in the {\it supervised} setting and attains the semi-parametric efficiency bound under a fully non-parametric model van2000asymptotic for estimating $\theta_0$. We will therefore use the limiting supervised posterior $\mathcal{N}(\widehat \theta_{\sup}, \sigma^2_{Y}/n)$ as a {\it benchmark} for asymptotic estimation/inference efficiency comparisons with BDMI later.

A motivating imputation-type Bayesian SS approach

The construction of the supervised posterior $\Pi_{\sup}$ (and $\widehat\theta_{\sup}$) naturally does not utilize the large unlabeled data $\mathcal{U}$ on $\mathbf{X}$ available in the SS setting. By virtue of its large size, $\mathcal{U}$ essentially informs us on the distribution, $\mathbb{P}_{\mathbf{X}}$, of $\mathbf{X}$. Thus, whenever $\mathbb{P}_{\mathbf{X}}$ is informative about the parameter of interest zhang2000value, seeger2000learning, one may hope to utilize $\mathcal{U}$ and come up with an improved SS Bayesian estimation procedure with a more efficient, i.e., tighter posterior contracting around $\theta_0$ (albeit at a $\sqrt{n}$-rate, since information on $Y$ is still limited to $n$ observations), and accordingly a $\sqrt{n}$-consistent point estimator of $\theta_0$ that is more efficient than $\widehat\theta_{\sup}$. We now discuss such an intuitive {\it imputation-based} approach with a natural Bayesian flavor, along with its potential drawbacks, which form a crucial basis for our final formulation of the BDMI method in Section (ref).

Recalling $m_0(\mathbf{X}) \equiv \mathbb{E}(Y \mid \mathbf{X})$, the functional $\theta_0 \equiv \theta_0(\mathbb{P}_{\mathbf{Z}}) = \mathbb{E}(Y)$ can be written via iterated expectations as: $\theta_0 \equiv \theta_0(\mathbb{P}_{\mathbf{X}}; m_0) = \mathbb{E}_{\mathbf{X}}\{\mathbb{E}_{Y \mid \mathbf{X}}(Y \mid \mathbf{X})\} =\mathbb{E}_{\mathbf{X}}\{m_0(\mathbf{X})\}$. This representation clearly explains the {\it connection} between $\mathbb{P}_{\mathbf{X}}$ and $\theta_0$, and the potential for $\mathcal{U}$ to be exploited through bringing in the {\it nuisance} function $m_0$ (unknown but {\it estimable} via $\mathcal{L}$). One can then construct an imputation-based Bayesian SS approach as follows.

Suppose one learns $m_0(\cdot)$ from $\mathcal{L}$ via {\it any} reasonable Bayesian regression method (see Remark (ref) for some examples) that provides a {\it nuisance posterior} $\Pi_{\mathbf{m}} \equiv \Pi_{\mathbf{m}}(\cdot\,;\mathcal{L})$ for $m$. Then, using the identity $\theta_0 = \mathbb{E}_{\mathbf{X}}\{m_0(\mathbf{X})\}$, and replacing $\mathbb{E}_{\mathbf{X}}$ therein with an empirical average over $\mathcal{U}$, one may obtain an {\it induced posterior} $\Pi_{\text{imp}}$ for $\theta$ via a natural {\it imputation} approach, i.e., for samples $\widetilde m \sim \Pi_{\mathbf{m}}$, we let $\theta_{\text{imp}} \equiv \theta_{\text{imp}}(\widetilde m) := {N^{-1}}\sum_{\mathbf{X}_i \in \mathcal{U}} \widetilde m (\mathbf{X}_i) \sim \Pi_{\text{imp}}$. Further, by linearity of expectation, it is easy to show that, $\widehat\theta_{\text{imp}}:= N^{-1} \sum_{\mathbf{X}_i \in \mathcal{U}} \widehat m(\mathbf{X}_i)$ is the posterior mean of $\Pi_{\text{imp}}$ (and hence, a point estimate of $\theta_0$), where $\widehat m(\cdot) :=\mathbb{E}_{\widetilde m \sim \Pi_{\mathbf{m}}}\{\widetilde m(\cdot) \mid \mathcal{L}\}$ is the posterior mean of $\Pi_{\mathbf{m}}$.

There are two major issues with this approach: (i) potential {\it misspecification} of $\Pi_{\mathbf{m}}$ in learning the true $m_0$; and (ii) more importantly, effect of the {\it nuisance} $\Pi_{\mathbf{m}}$'s {\it first-order properties (its rate/bias and variability) directly impacting} the {\it target} $\Pi_{\text{imp}}$'s {\it first-order behavior}. To illustrate, consider the {\it ideal} case: $N = \infty$. Then, the posterior sample $\theta_\text{imp}$ equals $\mathbb{E}_{\mathbf{X}}\{\widetilde m(\mathbf{X}) \hspace{0.01in} | \hspace{0.01in} \widetilde m\} \equiv \theta_0 + \mathbb{E}_{\mathbf{X}}\{\widetilde m(\mathbf{X}) - m_0(\mathbf{X}) \hspace{0.01in} | \hspace{0.01in} \widetilde m\}$. Thus, when misspecification is allowed, i.e., $\mathbb{E}_{\widetilde m \sim \Pi_{\mathbf{m}}} \{ \|\widetilde m(\mathbf{X}) - m^*(\mathbf{X}) \|_{\mathbb{L}_2(\mathbb{P}_{\mathbf{X}})} \mid \mathcal{L}\} \overset{\mathbb{P}}\to 0$ (under $\mathbb{P}_{\mathcal{L}}$) for {\it some} function $m^*(\cdot) \in \mathbb{L}_2(\mathbb{P}_{\mathbf{X}})$ possibly $\neq m_0(\cdot)$, then $\Pi_{\text{imp}}$ may become {\it inconsistent} (i.e., not contracting around the true $\theta_0$). More fundamentally, {\it even if} $m^*(\cdot) = m_0(\cdot)$, the {\it entire} first-order behavior (rate, shape, and variability) of $\Pi_{\text{imp}}$ depends {\it directly} on the corresponding behavior of (posterior of): $\widetilde m (\cdot) - m_0 (\cdot)$, the `bias term', making $\Pi_{\text{imp}}$ {\it sensitive}, in the {\it first order}, to $\Pi_{\mathbf{m}}$'s first order properties, and accordingly, the choice of the {\it method} used therein. In particular, if $\Pi_{\mathbf{m}}$ has a contraction rate, $a_n$, {\it slower} than $n^{-1/2}$, then so will $\Pi_{\text{imp}}$. More importantly, the {\it variability} of $\Pi_{\text{imp}}$ itself (after scaling by its rate) will be directly impacted by that of $\Pi_{\mathbf{m}}$. Overall, this indicates that to obtain a BvM-type result on $\Pi_{\text{imp}}$ -- necessary to ensure provably valid estimation and inference on $\theta_0$ -- one {\it requires} the availability of a corresponding semi-parametric BvM-type result under the nuisance $\Pi_{\mathbf{m}}$, which may necessitate delicate conditions/control on specifics of $\Pi_{\mathbf{m}}$'s construction. This becomes especially challenging when using non-smooth or complex methods, e.g., sparse regression in high dimensions or non-parametric machine learning methods, as nuisance estimators. These methods, while highly relevant and popular, have rates slower than $n^{-1/2}$, as well as unclear first-order properties with often intractable posteriors and limited availability (or feasibility) of corresponding BvM results. In general, this first-order sensitivity of $\Pi_{\text{imp}}$ and its reliance on such intricate aspects of $\Pi_{\mathbf{m}}$, therefore, jeopardizes rate-optimal and provably valid inference on $\theta_0$ with the correct variance. In Section (ref) of the \hyperref[sec:supplementary]{Supplementary Material}, we present a detailed case study on $\Pi_{\text{imp}}$ (and also compare it to BDMI) showcasing its sensitivity and failure to provide a valid inference on $\theta_0$.

Bayesian debiased modeling and inference: BDMI

This section introduces the BDMI approach, which addresses the limitations of the imputation approach discussed in Section (ref), by appropriately accounting for nuisance estimation bias within a Bayesian likelihood framework. BDMI is based on the principle of disentangling the nuisance parameter, and jointly {\it learning} its bias with the parameter of interest via targeted summary statistics {\it amenable} to Bayesian modeling. Incorporating this debiasing idea and the targeted modeling approach are our key methodological contributions towards {\it Bayesian semi-parametric inference}, in general, for robust and efficient inference in the presence of high dimensional nuisances, drawing parallels to the recent frequentist DML literature chernozhukov2018double.

Bayesian debiasing: Overcoming the bias from nuisance estimation within the Bayesian framework

For exposition of the BDMI approach and its salient features, we assume for the time being that there exists a dataset $\mathcal{S}$ which is an {\it independent copy} of the labeled data $\mathcal{L}$. The sample size $s_n$ of $\mathcal{S}$ is assumed to be of the same order as $n$; see Section (ref) for more details. Suppose the nuisance estimation is performed on this $\mathcal{S}$, using {\it any} reasonable Bayesian (or frequentist) method by constructing a likelihood for the nuisance parameter $m$ on $\mathcal{S}$, combining with a suitable prior on $m$, to obtain a posterior $\Pi_{\mathbf{m}}$ for $m$. For our primary goal of inference on $\theta_0$, the specific construction of $\Pi_{\mathbf{m}}$ is not crucial, provided it satisfies some basic regularity conditions (see Section (ref) for details). Henceforth, we assume access to a {\it generic} posterior $\Pi_{\mathbf{m}}$ for $m$, noting that $\Pi_{\mathbf{m}}(\cdot) \equiv \Pi_{\mathbf{m}}(\cdot; \mathcal{S})$ is itself a {\it random} distribution dependent on $\mathcal{S}$. For simplicity, this dependence is suppressed in the notations whenever clear from context. The dataset $\mathcal{S}$ can be viewed as {\it training data}, used {\it solely} to obtain the nuisance posterior $\Pi_{\mathbf{m}}$ for $m$. In contrast, $\mathcal{D} = \mathcal{L} \cup \mathcal{U}$ serves as {\it test data}, used to obtain the posterior for the parameter of interest $\theta$ via the BDMI procedure. In practice, we construct such pairs of independent training and test datasets from the original data $\mathcal{D}$ itself via {\it sample splitting}; see Section (ref).

Let $\widetilde{m} : \mathbb{R}^p \to \mathbb{R}$ be any {\it random function} van2000asymptotic output from $\mathcal{S}$ (e.g., a posterior sample from a Bayesian regression model fitted to $\mathcal{S}).$ More formally, $\widetilde{m}: (\Omega_{\mathcal{S}}, \mathbb{P}_{\mathcal{S}}) \times \mathbb{R}^p \to \mathbb{R}$ is a measurable map, i.e., $\{\widetilde{m}(x)\}_{x \in \mathbb{R}^p} \equiv \{\widetilde{m}(\omega;x)\}_{x \in \mathbb{R}^p}$ is a {\it stochastic process}, with sample paths $\widetilde{m} (\omega; \cdot)$ for $\omega \in \Omega_{\mathcal{S}}$, where $(\Omega_{\mathcal{S}}, \mathbb{P}_{\mathcal{S}})$ denotes the probability space underlying the randomness of $\mathcal{S}$ and any derived measures (e.g., posteriors) from it. Suppose now the argument $\mathbf{x}$ (or domain) of $\widetilde m(\mathbf{x}) \equiv \widetilde m(\omega; \mathbf{x})$ is {\it measurized} (randomized) {\it independently} as: $\mathbf{X} \sim \mathbb{P}_{\mathbf{X}} \perp \!\!\! \perp \mathbb{P}_{\mathcal{S}}$, e.g., $\mathbf{X} \in \mathcal{D}$ meets this requirement, since $\mathcal{D} \perp \!\!\! \perp \mathcal{S}$ by construction. Consider the {\it doubly} random variable $\widetilde{m}(\mathbf{X})$ -- having two sources of randomness that are independent -- (i) the process $\widetilde{m}(\cdot)$ {\it itself} from $\mathcal{S}$, and (ii) its random argument $\mathbf{X} \sim \mathbb{P}_{\mathbf{X}}$ from $\mathcal{D}$ ($\perp \!\!\! \perp \mathcal{S}$). We can then write $\theta_0 \equiv \mathbb{E}(Y)$ as:

The steps in both (ref)--(ref) use $\widetilde{m}(\cdot)$ from $\mathcal{S}$ is $\perp \!\!\! \perp$ of $\mathbf{X}$ (and $\mathbf{Z})$ $\in \mathcal{D}$. This independence is crucial and necessary to derive (ref), which we refer to as the {\it debiased} {\it representation} of $\theta_0$. For notational clarity, we emphasize that for a given $\widetilde{m}$, $b(\widetilde{m})$ should be interpreted as a parameter dependent on $\widetilde{m}$, i.e., a {\it function} of $\widetilde{m}$. Finally, we reiterate that the above representations (ref)--(ref) remain valid if $\widetilde{m}(\cdot)$ is a random draw from the posterior $\Pi_{\mathbf{m}}$, and $\mathbf{X} = \mathbf{X}_i \in \mathcal{D}$ ($i = 1, \ldots, n + N$), and $\mathbf{Z} = \mathbf{Z}_i \in \mathcal{L}$ ($i=1,\ldots,n$), since $\Pi_{\mathbf{m}}$ is constructed from $\mathcal{S}$ which is independent of $\mathcal{D}$. Subsequent references to (ref)--(ref) are with respect to (w.r.t.) these particular choices.

Note that the first term $b(\widetilde{m})$ in (ref) is essentially the expected {\it bias}, which is the price of replacing $m_0(\cdot)$ with a random sample $\widetilde{m}(\cdot)$. As noted in Section (ref), this is precisely the primary cause of the issues with the imputation approach. Modeling this $b(\widetilde{m})$ itself, along with $\theta_0$, is the central idea of BDMI. Note further that:

This shows $b(\widetilde{m})$ captures two pivotal aspects: (i) when $m^*(\cdot) \neq m_0(\cdot)$, the first term measures its average deviation from $m_0(\cdot)$, and (ii) the second term importantly reflects the {\it variability} of $\widetilde{m}(\cdot)$ itself as a sample from $\Pi_{\mathbf{m}}$ (which is further random through $\mathcal{S}$). From the perspective of statistical learning theory Vapnik1998, one could think of the first term as {\it approximation error} and the second term as {\it estimation error}.

Most importantly, observe that (ref) implies we also have {\it i.i.d. replicates} $\{Y_i - \widetilde{m}(\mathbf{X}_i)\}_{i \in \mathcal{L}}$ and $\{\widetilde{m}(\mathbf{X}_i)\}_{i \in \mathcal{U}}$ from {\it conditionally (given $\widetilde m$) independent sources} that {\it target} {$b(\widetilde{m})$} and {$\theta_0 - b(\widetilde m)$,} respectively, through their expectations. Thus, $b(\widetilde{m})$ and $\theta_0 - b(\widetilde m)$ can be seen as {\it functionals} of the underlying distribution of $\mathcal{L}$ and $\mathcal{U}$, specifically depending on the {\it summary statistics} (means) of $Y - \widetilde{m}(\mathbf{X})$ in $\mathcal{L}$ and $\widetilde{m}(\mathbf{X})$ in $\mathcal{U}$ (given $\widetilde{m}$ from an independent source), respectively. The {\it basic premise} of BDMI is: to model the data for these {\it target-specific parameters} -- $b(\widetilde{m})$ and $\theta_0 - b(\tilde m)$ -- via summary statistics, since they {\it directly inform us on $\theta_0$, while also learning the bias} induced by $\widetilde{m}$. This {\it targeted modeling} of summary statistics (instead of the entire data as in traditional Bayesian approaches) is a salient feature of BDMI. Further, its modeling of the bias $b(\widetilde m)$ {\it encodes a Bayesian form of debiasing} which plays a crucial role in ensuring nuisance-insensitive inference for $\theta_0$.

Targeted modeling of summary statistics: Likelihood construction and final posterior

We are now ready to introduce the target-specific model construction discussed in the previous section. Given $\widetilde{m} \sim \Pi_{\mathbf{m}}$ (from $\mathcal{S}$), the i.i.d. replicates $\{ Y_i - \widetilde{m}(\mathbf{X}_i) \}_{i = 1}^n$ and $\{\widetilde{m}(\mathbf{X}_i)\}_{i = n+1}^{n+N}$ from $\mathcal{D}$ ($\perp \!\!\! \perp \mathcal{S}$) target $b(\widetilde{m})$ and $\theta_0 - b(\widetilde{m})$, respectively, in terms of their means. These variables are now treated as our `observables' on the data $\mathcal{D} \mid \widetilde{m}$, and we now present a working likelihood construction for these observables on this data. To proceed, let us first define $\sigma^2_{1}(\widetilde{m}) := \mathrm{Var}_{\mathbf{Z}}\{Y - \widetilde{m}(\mathbf{X})\}$ and $\sigma^2_{2}(\widetilde{m}) := \mathrm{Var}_{\mathbf{X}}\{\widetilde{m}(\mathbf{X})\}$. Then, {\it given} $\widetilde{m}$, $Y_i - \widetilde{m}(\mathbf{X}_i)$ are i.i.d. with mean $b(\widetilde{m})$ and variance $\sigma^2_{1}(\widetilde{m})$ for $i \in \{1, \dots , n\}$, and $\widetilde{m}(\mathbf{X}_i)$ are i.i.d. with mean $\theta_0 - b(\widetilde{m})$ and variance $\sigma^2_{2}(\widetilde{m})$ for $i \in \{n + 1, \dots , n + N\}$. Since these observables are i.i.d., a natural choice of a working model for such data could be based on Normal distributions with unknown variances, as follows:

Then, the likelihood as a function of the parameters $\{\theta, b(\widetilde{m}), \sigma^2_{1}(\widetilde{m}), \sigma^2_{2}(\widetilde{m})\}$ is given by:

The (pseudo-) likelihood constructed above can be combined with a prior distribution on the model parameters $\{\theta, b(\widetilde{m}), \sigma^2_{1}(\widetilde{m}), \sigma^2_{2}(\widetilde{m})\}$ using Bayes' formula to yield a posterior, and thereafter a {\it marginal posterior} $\Pi_{\boldsymbol \theta}$ of $\theta$.

We note that the Normal distributions in (ref) above are only chosen as {\it working}, i.e., not necessarily correctly specified, distributions. Since a posterior depends on the data only through sufficient statistics, one could directly model the sample averages of $Y - \widetilde{m}(\mathbf{X})$ and $\widetilde{m}(\mathbf{X})$ as Normally distributed with appropriate parameters under modeling assumptions similar in spirit to (ref), operationally leading to the same posterior. In that case, one could simply treat the sample means as the `derived' observations, and since, given a sufficiently large number of observations, the sample averages are approximately Normal following the Central Limit Theorem (CLT), the Normality assumption on the sample averages would therefore be quite reasonable.

As a concrete {\it prior choice}, for the sake of theoretical and computational simplicity, we recommend using an improper prior on the model parameters $\{\theta, b(\widetilde{m})$, $\sigma^2_{1}(\widetilde{m}), \sigma^2_{2}(\widetilde{m})\}$ in (ref), given by:

with $\sigma^2_{1}(\widetilde{m})$ and $\sigma^2_{2}(\widetilde{m})$ being independent. We note that more general prior choices could also be employed here (see Remark (ref) for a discussion) without altering the asymptotic conclusions, such as the limiting posterior and related properties of the procedure, established in Section (ref). For instance, by defining $\delta := \theta - b(\widetilde{m})$, one could place independent conjugate Normal-Inverse Gamma priors on $\{b(\widetilde{m}), \sigma^2_{1}(\widetilde{m})\}$ and $\{\delta, \sigma^2_{2}(\widetilde{m})\}$. The proposed improper prior in (ref) can then be viewed as a limiting (diffused) version of such a proper prior.

We now explicitly compute the marginal posterior $\Pi_{\boldsymbol \theta}$ of $\theta$ under (ref) and the prior choice (ref), as follows.

Note that $\Pi_{\boldsymbol \theta}$, being a convolution of two $t$-distributions, is {\it easy to sample} from (e.g., for constructing CIs). Further, the {\it posterior mean:} $\widehat \theta_{\mathrm{BDM}}(\widetilde m)$ of $\Pi_{\boldsymbol \theta}$ can be considered as a natural point estimator of $\theta_0$. Note that the $\widetilde m$ in $\widehat \theta_{\mathrm{BDM}}(\widetilde m)$ reflects that the estimator (and the posterior $\Pi_{\boldsymbol \theta} \equiv \Pi_{\boldsymbol \theta}(\widetilde m)$ itself) fundamentally depends on the nuisance posterior sample $\widetilde m \sim \Pi_{\mathbf{m}}$ used. From Proposition (ref), it follows that $ \widehat \theta_{\mathrm{BDM}}(\widetilde m) = \mu_{n}(\widetilde{m}) + \mu_{N}(\widetilde{m})$. Note that $\widehat \theta_{\mathrm{BDM}}(\widetilde m)$ (and $\Pi_{\boldsymbol \theta}$, in general) utilize {\it both} $\mathcal{L}$ and $\mathcal{U}$, thereby justifying its billing as an SS approach. Also, as $n, N \to \infty$, it converges to $\theta_0$ {\it even if} $\Pi_{\mathbf{m}}$ is misspecified. This is because the first term in $\widehat \theta_{\mathrm{BDM}}(\widetilde m)$ targets $\mathbb{E}_{\mathbf{Z}} [\{Y - \widetilde{m}(\mathbf{X}) \}\mid \widetilde{m}]$, while the second term targets $\mathbb{E}_{\mathbf{X}}\{\widetilde{m}(\mathbf{X}) \mid \widetilde{m}\}$, hence canceling out $\widetilde m$'s effect. Thus, BDMI gives a posterior mean that is {\it always} a consistent point estimator. Moreover, one would expect the spread of the posterior $\Pi_{\boldsymbol \theta}$ to be of the correct rate $n^{-1/2}$, and also tighter than the supervised counterpart. These claims, along with other desirable properties of BDMI, are formally established later in Section (ref).

Sample splitting based version: BDMI with cross-fitting (BDMI-CF)

To practically implement the ideas introduced in Sections (ref) and (ref), we need to construct independent training and test dataset pairs $(\mathcal{S}, \mathcal{D})$ such that $\mathcal{S} \perp \!\!\! \perp \mathcal{D}$. To achieve this from the original data $\mathcal{D} = \mathcal{L} \cup \mathcal{U}$, we employ a $K$-fold sample splitting (with cross-fitting) procedure, where $K \geq 2$ is {\it fixed} (relative to $n,N$) and we assume without loss of generality (w.l.o.g.), that $|\mathcal{L}| = n$ and $|\mathcal{U}| = N$ are divisible by $K$. To construct independent training and test datasets required for the debiasing representation in (ref), we perform $K$-fold sample splitting by randomly partitioning the indices $\{1, \dots , n\}$ (for $\mathcal{L}$) and $\{ n+1, \dots n+N\}$ (for $\mathcal{U}$) into $K$ disjoint folds $\{\mathcal{I}_k\}_{i = 1}^K$ and $\{\mathcal{J}_k\}_{i = 1}^K$, respectively, with each fold $\mathcal{I}_k$ of size $n_K := n/K$ and $\mathcal{J}_k$ of size $N_K := N/K$, for each $k \in \{1, \dots , K\}$, define $\mathcal{I}_k^- := \{1, \dots , n\} \backslash I_k$. Then, using these partitions, we construct pairs of training and test data folds $\{(\mathcal{S}_k, \mathcal{D}_k)\}_{k = 1}^K$, where $\mathcal{S}_k := \{\mathbf{Z}_i: i \in \mathcal{I}_k^- \}$ $\perp \!\!\! \perp$ $\mathcal{D}_k := \mathcal{L}_k \cup \mathcal{U}_k$, with $\mathcal{L}_k := \{\mathbf{Z}_i: i \in \mathcal{I}_k \}$ and $\mathcal{U}_k := \{\mathbf{X}_i: i \in \mathcal{J}_k \}$. This provides $K$ such (training, test) data pairs for constructing the BDMI approach on each pair. Importantly, the test datasets $\mathcal{D}_1, \dots, \mathcal{D}_K$ are all disjoint and {\it independent}.

Adopting the BDMI construction from Section (ref), we now detail the BDMI procedure for one pair $(\mathcal{S}_k, \mathcal{D}_k)$. Since $\mathcal{S}_k \perp \!\!\! \perp \mathcal{D}_k$, we use the training subfold $\mathcal{S}_k$ to obtain the nuisance posterior $\Pi_{\mathbf{m}}^{(k)}$ for $m$, as detailed in Section (ref). Let $\widetilde{m}_k$ be {\it one} random sample from $\Pi_{\mathbf{m}}^{(k)}$. Following the same model construction in Section (ref), we use the same likelihood formulation for the test subfold $\mathcal{D}_k$ as given in equations (ref)--(ref):

Using the same improper prior on the model parameters $\{\theta, b(\widetilde{m}_k), \sigma^2_{1}(\widetilde{m}_k), \sigma^2_{2}(\widetilde{m}_k)\}$ from (ref), and applying Proposition (ref) with $(\mathcal{S},\mathcal{D})$ therein set as $ (\mathcal{S}_k, \mathcal{D}_k)$, we derive the marginal posterior $\Pi_{\boldsymbol \theta}^{(k)}$ for $\theta$ as follows:

Consistent with our earlier notation, let $\widehat \theta_{\mathrm{BDM}}^{(k)}(\widetilde{m}_k)$ denote the posterior mean of $\Pi_{\boldsymbol \theta}^{(k)}$. From Proposition (ref), we have $ \widehat \theta_{\mathrm{BDM}}^{(k)}(\widetilde{m}_k) = \mu_{n_K}(\widetilde{m}_k) + \mu_{N_K}(\widetilde{m}_k) $ and it retains the same properties as $\widehat \theta_{\mathrm{BDM}}(\widetilde{m})$ from Section (ref).

While sample splitting enables us to obtain the debiased representation in (ref), which is crucial for the BDMI approach, it uses only a subset $\mathcal{D}_k$ of the full dataset $\mathcal{D}$ to obtain a posterior for $\theta$. This causes a notable lack of efficiency. Since sample splitting produces $K$ splits, each data fold pair $(\mathcal{S}_k, \mathcal{D}_k)$ can be utilized to obtain a posterior $\Pi_{\boldsymbol \theta}^{(k)}$ of $\theta$ for $k = 1, \dots, K$. We now introduce a method for combining these posteriors of $\theta$, referred to as {\it BDMI with cross-fitting} (BDMI-CF), to construct an {\it aggregated} full-data posterior for $\theta$. This approach addresses the efficiency loss discussed earlier by fully utilizing the available data and ensuring that the variance and contraction rates of the final procedure depend directly on $n$, as shown in Theorem (ref).

BDMI-CF is inspired by the frequentist cross-fitting (CF) idea chernozhukov2018double, addressing challenges in high dimensional nuisance parameter estimation. The conventional CF approach has been used to (i) relax strong assumptions, e.g., Donsker class conditions van2000asymptotic, and (ii) make the sample splitting process efficient utilizing the full data in a `cross-fitted' manner chernozhukov2018double. CF techniques are widely used in the modern semi-parametric inference literature, where a combined estimator is obtained by averaging the estimators obtained from each split to regain full efficiency chernozhukov2018double, newey2018cross. In a {\it Bayesian} framework, however, additional care is required during the combination step, since entire {\it distributions (posteriors)} must be aggregated rather than point estimates. BDMI-CF addresses this issue by employing a consensus Monte Carlo-type approach Scott02042016 to suitably aggregate the posteriors from the sub-folds. This type of usage of cross-fitting (CF) for combining posteriors in Bayesian semi-parametric inference problems is not common. In the existing Bayesian literature, sample splitting has primarily been used to improve computational efficiency when handling large datasets scott2022bayes. However, BDMI leverages sample splitting in a novel way: to ensure independence between the estimation of the nuisance parameter and the parameter of interest, and further via CF based aggregation, ensures efficient usage of the {\it entire} data. We now discuss the CF procedure.

Let $\theta_1, \dots , \theta_K $ be independent random variables drawn from the corresponding posteriors $\Pi_{\boldsymbol \theta}^{(1)}, \dots, \Pi_{\boldsymbol \theta}^{(K)}$ which are obtained from $(\mathcal{S}_1,\mathcal{D}_1), \dots , (\mathcal{S}_K,\mathcal{D}_K)$, respectively. We then define a new random variable:

The distribution $\Pi_{\boldsymbol \theta}$ in (ref) is referred to as the {\it final (aggregated)} posterior of $\theta$ from BDMI, specifically BDMI-CF. This final posterior $\Pi_{\boldsymbol \theta}$ is a (scaled) {\it convolution} of the posteriors $\Pi_{\boldsymbol \theta}^{(1)}, \dots, \Pi_{\boldsymbol \theta}^{(K)}$ obtained from each data fold pair $(\mathcal{S}_1,\mathcal{D}_1), \dots , (\mathcal{S}_K,\mathcal{D}_K)$. Hence, samples from $\Pi_{\boldsymbol \theta}$ can be easily generated by construction.

Further, by linearity of expectation, the posterior mean $\widehat \theta_{\mathrm{BDM}}(\widetilde{m}_{\mathrm{CF}})$ of $\Pi_{\boldsymbol \theta}$ is the average of the posterior means $\widehat \theta_{\mathrm{BDM}}^{(1)}(\widetilde{m}_1), \dots, \widehat \theta_{\mathrm{BDM}}^{(K)}(\widetilde{m}_K)$ from the corresponding posteriors $\Pi_{\boldsymbol \theta}^{(1)}, \dots, \Pi_{\boldsymbol \theta}^{(K)}$. More explicitly,

where $\widetilde{m}_{\mathrm{CF}}(\mathbf{X}_i) := \widetilde{m}_k(\mathbf{X}_i)$ for $i \in \mathcal{I}_k$ or $i \in \mathcal{J}_k$ where $\widetilde{m}_k$ is a random sample from the respective posterior $\Pi_{\mathbf{m}}^{(k)}$ of $m$ for $k = 1, \dots, K$. Naturally, we consider the posterior mean $\widehat \theta_{\mathrm{BDM}}(\widetilde{m}_{\mathrm{CF}})$ as a point estimator of $\theta_0$. Furthermore, Theorem (ref) guarantees the $\sqrt{n}$-consistency of $\widehat \theta_{\mathrm{BDM}}(\widetilde{m}_{\mathrm{CF}})$ as an estimator of $\theta_0$. Detailed properties of $\widehat \theta_{\mathrm{BDM}}(\widetilde{m}_{\mathrm{CF}})$, and more generally the posterior $\Pi_{\boldsymbol \theta}$ in (ref), are further examined in Section (ref). We now present the final algorithm for our BDMI (specifically, BDMI-CF) approach in Algorithm (ref).

Theoretical properties of the BDMI procedure

In this section, we analyze in detail the theoretical underpinnings of our proposed BDMI procedure. Under mild regularity conditions, we show (in Theorems (ref)--(ref)) that the BDMI posteriors $\{\Pi_{\boldsymbol \theta}^{(k)}\}_{k = 1}^K$ (the `one fold' versions) and $\Pi_{\boldsymbol \theta}$ (the final aggregated version via CF) all inherit BvM-type limiting behaviors with asymptotically Gaussian posteriors contracting around the true $\theta_0$ at a $\sqrt{n}$-rate, along with various desirable properties on robustness, efficiency and nuisance insensitivity, which are all discussed in detail subsequently.

The BDMI-CF procedure provides the posterior $\Pi_{\boldsymbol \theta}$ with the posterior mean $\widehat\theta_{\mathrm{BDM}}(\widetilde{m}_{\mathrm{CF}})$ as defined in (ref). Naturally, $\widehat\theta_{\mathrm{BDM}}(\widetilde{m}_{\mathrm{CF}})$ can be considered as a valid {\it SS point estimator} for $\theta_0$. Beyond direct implications of Theorem (ref), the asymptotic behavior of the SS estimator $\widehat\theta_{\mathrm{BDM}}(\widetilde{m}_{\mathrm{CF}})$ inherently is of separate interest. Towards that, in Corollary (ref), we rigorously establish an asymptotically linear representation of $\widehat\theta_{\mathrm{BDM}}(\widetilde{m}_{\mathrm{CF}})$.

Robustness, efficiency and nuisance insensitivity of BDMI

Theorem (ref) establishes that, under the SS setting, the posterior $\Pi_{\boldsymbol \theta}$ concentrates around the true parameter $\theta_0$ at the parametric rate $1/ \sqrt{n}$ (ensuring usage of the {\it full} data) and possesses {\it universal robustness} to the choice of the nuisance estimation method. This robustness manifests in two ways: (i) {\it global robustness} w.r.t. the limiting function $m^*(\cdot)$, ensuring that $\Pi_{\boldsymbol \theta}$ contracts around $\theta_0$ at a rate $1/\sqrt{n}$ {\it regardless} of the contraction rate $a_n$ of $\Pi_{\mathbf{m}}$ and {\it even if} $m^*(\cdot) \neq m_0(\cdot)$; and (ii) {\it insensitivity} to the nuisance estimation bias, as $\Pi_{\boldsymbol \theta}$ is {\it not} affected by slower convergence rates $a_n$ of $\Pi_{\mathbf{m}}$, nor by $\Pi_{\mathbf{m}}$'s {\it own} first order properties like its shape, variability etc. (even after scaling by $a_n$). $\Pi_{\boldsymbol \theta}$ depends on $\Pi_{\mathbf{m}}$ {\it only} through its limit $m^*$, and validity/properties of $\Pi_{\boldsymbol \theta}$ as in Theorem (ref) requires only $a_n \to 0$. Hence, BDMI effectively addresses the primary issue of the imputation approach (see Section (ref)), where nuisance estimation bias directly characterizes the first-order behavior/properties of the posterior for $\theta$, and offers substantial {\it flexibility} in choosing regression methods to obtain $\Pi_{\mathbf{m}}$. In particular, it paves the way for using non-smooth or complex methods, like sparse regression (in high dimensions) or non-parametric ML methods, both of which may unavoidably have slow or unclear first order behaviors (refer to Remarks (ref) and (ref) for examples of these methods and their contraction rates). Moreover, BDMI-CF achieves {\it efficiency improvement} over the supervised approach based on $\mathcal{L}$, irrespective of whether $m^*(\cdot) = m_0(\cdot)$. While both $\Pi_{\boldsymbol \theta}$ and $\Pi_{\sup}$ converge to $\theta_0$ at the parametric rate $1/ \sqrt{n}$, the variance $\tau^2_{n, N}(m^*)$ of the limiting distribution is {\it always smaller} than the variance of the supervised approach as we will show in Remark (ref), and further achieves the {\it semi-parametric efficiency bound} when $m^*(\cdot) = m_0(\cdot)$ (correctly specified case). These results align with frequentist asymptotic theory in recent SS inference literature zhang2019semi, zhang2022high. Moreover, these desirable properties of $\Pi_{\boldsymbol \theta}$ also naturally extend to posterior summaries. In particular, the {\it posterior mean} $\widehat\theta_{\mathrm{BDM}}(\widetilde{m}_\mathrm{CF})$, as a valid SS point estimator of $\theta_0$, inherits these properties. As Corollary (ref) shows, it remains $\sqrt{n}$-consistent, asymptotically Normal, and asymptotically linear regardless of the nuisance estimation method, and its expansion is unaffected by the estimation bias/error of the nuisance, showing its first-order insensitivity. Finally, its asymptotic variance also equals the posterior variance $\tau^2_{n, N}(m^*)$ (see Remark (ref) below), ensuring {\it valid and accurate inference} for $\theta_0$.

A hierarchical variant of BDMI: h-BDMI

Recall that the original BDMI procedure, as described in Section (ref), is constructed using a {\it single} random sample $\widetilde{m} \sim \Pi_{\mathbf{m}}$. Alternatively, a more traditional Bayesian approach can be adapted by considering multiple samples of $\widetilde{m}$ through a hierarchical construction, as briefly mentioned in Remark (ref). This section presents this alternative version of BDMI, referred to as the {\it hierarchical-BDMI} (henceforth h-BDMI), which constructs a joint posterior of $(\theta, m)$ and then marginalizes over $m$ to obtain the marginal posterior of $\theta$. This differs from the original BDMI procedure, and h-BDMI aligns more closely with traditional hierarchical Bayesian modeling principles. For its exposition, we focus on only one data fold, say $\widetilde\mathcal{D}_k := \mathcal{D}_k \cup S_k$, where $\mathcal{D}_k$ and $\mathcal{S}_k$ are as defined in Section (ref) for some $k = 1, \dots, K$. Following the conventional Bayesian idea of integrating out the nuisance parameter $m$, we proceed as follows. Using $\mathcal{S}_k$ as a training data, we obtain a posterior $\Pi_\mathbf{m}^{(k)} \equiv \Pi_\mathbf{m}^{(k)}(\cdot; \mathcal{S}_k) $ for $m$. By the conditional independence between $m \sim \Pi_\mathbf{m}^{(k)}$ and $\mathcal{D}_k$, the joint posterior of $(\theta, m)$ has the pdf $ \pi(\theta, m \mid \widetilde\mathcal{D}_k) \ = \ \pi(\theta \mid m, \mathcal{D}_k) \hspace{0.7mm}\pi_{\mathbf{m}}^{(k)}(m), $ where $\pi_{\mathbf{m}}^{(k)}(\cdot)$ is the pdf of the nuisance posterior $\Pi_\mathbf{m}^{(k)}$ of $m$. The pdfs $\pi(\theta \mid m, \mathcal{D}_k)$ and $\pi_{\mathbf{m}}^{(k)}(m)$ remain as defined in Section (ref). By integrating out $m$, we obtain the marginal posterior of $\theta$, denoted $\widetilde \Pi_{\boldsymbol \theta}^{(k)}$, with corresponding pdf $\pi_{\boldsymbol \theta}^{(k)}(\cdot)$, based on h-BDMI as follows: $$ \pi_{\boldsymbol \theta}^{(k)}(\theta) \ = \ \int \pi\big(\theta, m \mid \widetilde\mathcal{D}_k\big) \hspace{0.4mm}\mathrm{d} m \ = \ \int \pi\big(\theta \mid m, \mathcal{D}_k\big)\hspace{0.4mm} \pi_{\mathbf{m}}^{(k)}(m) \hspace{0.4mm} \mathrm{d} m. $$ Estimation and inference on the true parameter $\theta_0 \equiv \mathbb{E}(Y)$ using h-BDMI can be performed based on this posterior $\widetilde \Pi_{\boldsymbol \theta}^{(k)}$, for any $k = 1,\ldots, K$. Using iterated expectations, the posterior mean $\widehat \theta_{\mathrm{hBDM}}^{(k)} \equiv \mathbb{E}_{\theta \sim \widetilde \Pi_{\boldsymbol \theta}^{(k)}}(\theta)$ of $\widetilde \Pi_{\boldsymbol \theta}^{(k)}$ can be expressed as $ \widehat \theta_{\mathrm{hBDM}}^{(k)} ~\equiv ~ \int \left\{ \int \theta \hspace{0.5mm} \pi(\theta \mid m, \mathcal{D}_k) \hspace{0.5mm} \mathrm{d} \theta \right\} \pi_{\mathbf{m}}^{(k)} \hspace{0.5mm} \mathrm{d} m, $ where the inner integral is the conditional mean of $\theta$ given $m$ and $\mathcal{D}_k$, i.e., $\mathbb{E}(\theta \mid m, \mathcal{D}_k)$. Under the prior choice in (ref), $\widehat \theta_{\mathrm{hBDM}}^{(k)}$ is explicitly given by: $$ \widehat \theta_{\mathrm{hBDM}}^{(k)} ~\equiv ~ \frac{1}{n_K} \sum_{i \in \mathcal{I}_k} \big\{Y_i - \widehat m^{(k)}(\mathbf{X}_i)\big\} + \frac{1}{N_K} \sum_{i \in \mathcal{J}_k} \widehat m^{(k)}(\mathbf{X}_i), ~\text{ where $\widehat m^{(k)}$ is the posterior mean of $\Pi_{\mathbf{m}}^{(k)}$}. $$ Also, it is easy to draw samples from the posterior $\widetilde \Pi_{\boldsymbol \theta}^{(k)}$ of $\theta$ to construct credible intervals. Specifically, for sufficiently large $M$, we first draw samples $\widetilde{m}_1, \dots , \widetilde{m}_M$ from the posterior $\Pi_{\mathbf{m}}^{(k)}$, and for each sample, we draw a sample $\theta \mid \widetilde{m}_j, \mathcal{D}_k$ from the posterior $\Pi^{(k)}_{\boldsymbol \theta}$ as described in Proposition (ref) for $j = 1, \dots , M$. This process yields $M$ samples of $\theta$ from the posterior $\widetilde \Pi_{\boldsymbol \theta}^{(k)}$. Finally, applying the h-BDMI procedure to each $\widetilde \mathcal{D}_1 , \dots , \widetilde \mathcal{D}_K$, we obtain the corresponding posteriors $\widetilde \Pi_{\boldsymbol \theta}^{(1)}, \dots \widetilde \Pi_{\boldsymbol \theta}^{(K)}$. Following the aggregation approach detailed in Section (ref), we can construct a CF-based aggregated posterior $\widetilde \Pi_{\boldsymbol \theta}$. These modifications can be incorporated into Algorithm (ref), which we omit for brevity. We next present the result on the theoretical properties of h-BDMI on one data fold $\widetilde\mathcal{D}_k$, followed by a discussion on the differences between BDMI and h-BDMI.

We conclude this section with a brief comparison between the BDMI and h-BDMI approaches. Firstly, Theorem (ref) establishes a corresponding BvM-type result for h-BDMI, similar to Theorem (ref) for the `one data fold' version of the original BDMI procedure described in Section (ref). While both theorems demonstrate that the marginal posteriors of $\theta$ inherit a BvM-type limiting behavior with the same limiting posterior, they {\it do} have some important differences. Notably, Theorem (ref) requires a {\it stronger} $L_1$-type (Bayes risk) convergence condition on the contraction of the posterior $\widetilde \Pi_{\mathbf{m}}^{(k)}$ around $m^*(\cdot)$, while Theorem (ref) relies on the much weaker in-probability type condition (ref). In practice, our simulation results in Section (ref) reveal that the difference between BDMI and h-BDMI is less pronounced. In most cases, the two methods perform similarly in estimating $\theta_0$, as illustrated in Table (ref). Occasionally, h-BDMI tends to give slightly conservative coverages compared to BDMI (see Table (ref)), which is not unexpected since h-BDMI involves multiple samples (hence more noise) as it integrates out the nuisance parameter $m$ rather than conditioning on a single draw.

A key advantage of BDMI lies in its simplicity and computational efficiency. Unlike h-BDMI, which requires multiple samples from the nuisance posterior $\Pi_{\mathbf{m}}$ of $m$, BDMI relies on only a single sample, reducing computation burden. Thus, we recommend the original BDMI approach for achieving {\it both} efficient estimation and reliable inference for the true parameter $\theta_0$. For further details and discussions, we refer to Section (ref).

Finally, while we have used a `one fold' version of h-BDMI here for clarity, it also admits a CF-based full data version (`h-BDMI-CF', if we may) analogous to the BDMI-CF procedure in Section (ref). This version inherits similar theoretical properties as Theorem (ref) (with the same distinctions as above). In our simulations in Section (ref), we implemented h-BDMI via this CF-based full data version to ensure a fair comparison with the BDMI-CF and supervised approaches. The notation `h-BDMI' therein refers to this CF-based version.

Numerical studies

We conducted extensive simulation studies to investigate the finite sample performance, both in estimation and inference, for our proposed SS approach(es) and the supervised approach under various settings. In particular, as point estimators, we compare the supervised estimator $\widehat\theta_{\sup} \equiv \overline{Y}$ based on $\mathcal{L}$, the posterior mean $\widehat\theta_{\mathrm{BDM}} \equiv \widehat\theta_{\mathrm{BDM}}(\widetilde{m}_{\mathrm{CF}})$ of $\Pi_{\boldsymbol \theta}$ from the final BDMI-CF procedure (as in Algorithm (ref)) and the posterior mean $\widehat\theta_{\mathrm{hBDM}}$ of $\widetilde \Pi_{\boldsymbol \theta}$ from the h-BDMI procedure (its CF based version) discussed in Section (ref). We compare their estimation efficiencies based on the empirical mean squared error ({\bf MSE}) and report their relative efficiencies ({\bf RE}) compared to the supervised estimator $\widehat\theta_{\sup}$. Further, for evaluating the accuracy of inference, we report the empirical coverage probabilities (CovP) and lengths ({\bf Len}) of the $95\%$ credible intervals ({\bf CI}s) obtained from their respective posteriors. Finally, as a performance benchmark for estimation efficiency, we also report the maximum (oracle) asymptotic relative efficiency ({\bf ORE}) relative to $\widehat\theta_{\sup}$, given by $\mathrm{Var}(\widehat\theta_{\sup})/\tau^2_{n, N}(m^*)$, where $\tau^2_{n, N}(m^*) = \mathrm{Var}\{Y - m^*(\mathbf{X})\}/n + \mathrm{Var}\{m^*(\mathbf{X})\}/N$ with $m^*(\cdot) = m_0(\cdot)$. For the choice of the number of folds $K$, we consider $K = 5$ and $10$. The reported simulation results are all based on 500 replications. We examine various true data generating mechanisms and different methods for nuisance parameter estimation, leading to both correctly specified and misspecified models for $m_0(\cdot)$. We discuss the correctly specified and misspecified model settings and their corresponding results in Sections (ref)-- (ref).

Simulation studies: Correctly specified models

Throughout, we set $n = 500$ and $N = 10000$, and considered $p = 50$ and $p = 166$ ($\approx n/3)$, representing moderate and high dimensional settings (relative to $n$), respectively. We generated $\mathbf{X} \sim \mathcal{N}_p(\mathbf{0}_p, I_p)$, and given $\mathbf{X} = \mathbf{x}$, we generated $Y \sim \mathcal{N}(m_0(\mathbf{x}), \sigma_0^2)$, where $m_0(\mathbf{x}) = \alpha_0 + \mathbf{x}'\boldsymbol\beta_0$ and $\sigma_0^2 = \mathrm{Var}\{m_0(\mathbf{X})\}/5$, and we used $\alpha_0 = 5$ and $\boldsymbol\beta_0 = (\mathbf{1}_{s/2}',\ \mathbf{0.5}_{s/2}',\ \mathbf{0}_{p-s}')'$ (for different choices of $s$ discussed below). Here, $\mathcal{N}_d(\boldsymbol{\mu}, \boldsymbol{\Sigma}$) denotes the $d$-variate ($d \geq 2$) Gaussian distribution with mean $\boldsymbol{\mu}_{d \times 1}$ and covariance matrix $\boldsymbol{\Sigma}_{d \times d}$, $I_d$ denotes the identity matrix of order $d$, and the notation $\mathbf{a}_{l}$, for any positive integer $l$ (e.g., $l = p$, $s/2$ or $p - s$, as above), denotes the vector $(a, \dots, a)'_{l \times 1}$ for any $a \in \mathbb{R}$ (e.g., $a = 0$, $0.5$ or $1$, as above). The parameter $s$ in $\boldsymbol\beta_0$ above denotes the {\it sparsity} of $\boldsymbol\beta_0$. For $p = 50$, we set $s = 7$ $(\approx \sqrt{p})$, or $s = 50 \equiv p$; while for $p = 166$, we set $s = 13$ ($\approx \sqrt{p}$), $ s= 55$ ($\approx p/3$), $s = 83$ ($\approx p/2$), or $s = 166 \equiv p$. These choices of $s$ span a variety of settings, including {\it sparse} $(s= \sqrt{p}$), {\it moderately dense} $(s = p/2$ or $p/3$), or {\it fully dense} ($s =p$) cases. Note that, except for the sparse case, none of these choices correspond to settings where $s$ (or $p$) may be considered small or fixed relative to $n$, and therefore appropriate sparsity-friendly nuisance estimation methods may {\it still} fail to consistently estimate $m_0$. For illustrative purposes, we consider {\it three choices} (all parametric model based) for obtaining the nuisance posterior $\Pi_{\mathbf{m}}$: Bayesian ordinary linear regression ({\tt Bols}), Bayesian ridge regression ({\tt Bridge}), and a sparse Bayesian linear regression method ({\tt Bsparse}) based on non-local priors (NLP) johnson2012bayesian. In all cases, we consider the Gaussian linear regression model $Y_i \mid \mathbf{X}_i, \alpha, \boldsymbol\beta, \sigma \overset{\text{i.i.d.}}\sim \mathcal{N}(\alpha +\mathbf{X}_i'\boldsymbol\beta, \sigma^2)$ for $i = 1, \dots, n$. For {\tt Bols}, we use a prior on $(\alpha,\boldsymbol\beta, \sigma^2)$ given by: $\pi(\alpha,\boldsymbol\beta \mid \sigma^2) \propto 1$ and $\pi(\sigma^2) \propto (\sigma^2)^{-1}$; and for {\tt Bridge}, the prior employed on $(\alpha, \boldsymbol\beta,\sigma^2)$ is: $\pi(\alpha \mid \sigma^2) \propto 1$, $\boldsymbol\beta \mid \lambda, \sigma^2 \sim \mathcal{N}_p(\mathbf{0}_p, \lambda^{-1}\sigma^2 I_p)$, with $\alpha$ and $\beta$ being independent, and $\pi(\sigma^2) \propto (\sigma^2)^{-1}$. We use an empirical Bayes approach to plug in a point estimate $\widehat\lambda$ for the prior precision (or ridge) parameter $\lambda$. The estimate $\widehat\lambda$ is obtained from the {\tt R} package {\tt glmnet} so that the posterior mean of $(\alpha, \boldsymbol\beta')' \in \mathbb{R}^{(p+1)}$ coincides with the cross-validated point estimate obtained from {\tt cv.glmnet} in the {\tt glmnet} package. For both these methods, we obtain that the posteriors of $(\alpha, \boldsymbol\beta)$ are multivariate $t$-distributions. For the {\tt Bsparse} method, we use the R package mombf to obtain posterior samples for $(\alpha, \boldsymbol\beta)$. The implementation details of the mombf and {\tt glmnet} packages are provided in Section (ref) of the \hyperref[sec:supplementary]{Supplementary Material}.

Table (ref) and Tables (ref)--(ref) present the results on estimation efficiency and inference, respectively, along with illustrations of the posteriors and their overall behaviors in Figures (ref)--(ref). As seen from Table (ref) (as well as the box plots in Figures (ref)--(ref)), the REs of $\widehat\theta_{\mathrm{BDM}}$ and $\widehat\theta_{\mathrm{hBDM}}$ w.r.t. $\widehat\theta_{\sup}$, i.e., $\text{MSE}(\widehat\theta_{\sup})/ \text{MSE}(\widehat\theta_{\mathrm{BDM}})$ and $\text{MSE}(\widehat\theta_{\sup})/ \text{MSE}(\widehat\theta_{\mathrm{hBDM}})$, are consistently greater than 1, ranging roughly between 2 to 5 across most settings. This highlights the substantial efficiency improvement achieved by BDMI over the supervised approach. In addition, as illustrated in Figures (ref)--(ref), apart from point estimators, the {\it posteriors themselves} are consistently and significantly {\it tighter} than the supervised posteriors, while throughout resembling a Gaussian behavior centered at the true $\theta_0$. These patterns hold generally {\it regardless} of the setting and/or the nuisance posterior.

Furthermore, Table (ref) illustrates that the efficiency improvement depends primarily on the dimensionality $p$ and the sparsity level $s$. In moderate-dimensional settings ($p = 50$), BDMI achieves (near-)optimal efficiency gains, with RE values close to each other and approaching the ORE value, regardless of the sparsity levels (sparse $s = \sqrt{p}$ and fully dense $s=p$). This confirms that BDMI performs {\it optimally} when the model is correctly specified {\it and} estimated well enough. The impact of the sparsity level becomes particularly apparent in high-dimensional scenarios ($n = 500$ with $p = 166$), where finite-sample nuisance estimation {\it bias} introduces a {\it soft} form of misspecification. Specifically, sparsity-friendly nuisance models (e.g., Bsparse) struggle to consistently estimate $m_0(\cdot)$ in moderately dense ($s = p/2$) or fully dense ($s = p$) settings, leading to somewhat lower RE values. However, in sparse settings $s = \sqrt{p}$ ($s = 13$), Bsparse achieves RE values that are close to ORE by leveraging the underlying sparse structure, outperforming non-sparse methods such as Bols and Bridge. Conversely, in fully dense settings ($p = s$, $s = 166$), Bsparse struggles to adapt and estimates a nearly constant function, resulting in RE values close to 1. In contrast, the non-sparse methods Bols and \texttt{Bridge} still target non-trivial approximations of $m_0$, yielding reasonably high RE values (approximately 3.70; see Table (ref)). These observations highlight that while BDMI remains robust under soft misspecification, the choice of a nuisance model can influence the {\it extent} of efficiency gain in dense settings, emphasizing the interplay among both $p$ and $s$. Notably, even in high-dimensional (fully dense) cases, RE values remain still acceptable (RE $> 1$), albeit not optimal, and BDMI consistently provides correct coverage around $95\%$ regardless of the nuisance parameter estimation methods. Moreover, Table (ref) shows that the RE values of the SS estimators tend to be slightly higher for $K = 10$. But in general, the results -- both for estimation and inference -- seem to be fairly robust across both choices of $K$. We thus recommend either choice in practice.

Tables (ref)--(ref) exhibit that BDMI {\it consistently} achieves {\it correct} coverage probabilities for $\theta_0$, maintaining approximately $95\%$ coverage across all settings with various choices of $p, s, K$, as well as different methods for obtaining the nuisance posterior $\Pi_{\mathbf{m}}$. This highlights the {\it robustness} of BDMI in providing valid and accurate inference (correct coverage), as well as substantial {\it improvement} over supervised inference with {\it tighter} CIs (typically around 50% tighter) across settings -- thereby validating its construction and our claimed theoretical properties. Figures (ref)--(ref) provide visual confirmation of these findings, showing that BDMI-based posteriors consistently exhibit always tighter spread than the supervised posterior, regardless of the setting or the nuisance posterior method. Additionally, the variability of the BDMI posteriors {\it remains} consistent within each setting, further emphasizing its robustness and nuisance-insensitivity across the different scenarios.

Lastly, comparing the BDMI and h-BDMI approaches, we observe that despite the former requiring only one sample from $\Pi_{\mathbf{m}}$, both methods perform similarly across most settings, which: (i) validates our earlier claims on their common theoretical properties, {and} (ii) also {\it reinforces} the crucial role of {\it debiasing} common to both, that {\it negates} any distinction between the use of one vs. many $\widetilde m$ samples. The point estimators $\widehat\theta_{\mathrm{BDM}}$ and $\widehat\theta_{\mathrm{hBDM}}$ show very similar efficiencies with h-BDMI marginally higher in some cases, while for inference, h-BDMI often tends to give slightly conservative coverages $> 95\%$ (likely due to more noise from its hierarchical nature). Overall, given its computational simplicity, we recommend the original BDMI approach.

Simulation studies: Misspecified models

Section (ref) considered scenarios where the true model is linear, with Bayesian linear methods used to obtain the nuisance posterior $\Pi_{\mathbf{m}}$ of $m$. Although the models were technically “correctly" specified, high dimensional (and dense) settings do {\it not} necessarily guarantee consistent estimation of the true $m_0(\cdot)$, leading to a `soft' form of misspecification. We now examine the functional form of misspecification where the limiting function $m^*(\cdot)$ around which $\Pi_{\mathbf{m}}$ contracts, is {\it not} equal to $m_0(\cdot)$, i.e., $m_0(\cdot)$ is nonlinear but the {\it fitted} model for learning $\Pi_{\mathbf{m}}$ remains linear. Even in such cases, Theorems (ref)--(ref) ensure BDMI's validity with efficiency improvement persisting (see Table (ref)), though the improvements may not reach the optimal ORE.

Throughout, we set $N = 10000$ and $n = 500$. To illustrate our points, we specifically study non-linear, but low or moderate dimensional models with $p = 10$ (and sparsity $s = 10$ or $3$) or $p = 50$ (and sparsity $s = 50$ or $7$). We generated $\mathbf{X} \sim \mathcal{N}_p(\mathbf{0}, I_p)$ as in Section (ref) and given $\mathbf{X} = \mathbf{x}$, we generate $Y \sim \mathcal{N}(m_0(\mathbf{x}), \sigma_0^2)$ with $m_0(\mathbf{x}) = \alpha_0 + \mathbf{x}'\boldsymbol\beta_0 + (\mathbf{x}'\boldsymbol\gamma_0)^2$ and $\sigma_0^2 = \mathrm{Var}\{m_0(\mathbf{X})\}/5$. Here, $\alpha_0 = 5$, $\boldsymbol\beta_0 = (\mathbf{1}_{s/2}',\ \mathbf{0.5}_{s/2}',\ \mathbf{0}_{p-s}')'$ and $\boldsymbol\gamma_0$ is constructed to ensure $\sqrt{\mathbb{E}\{(\boldsymbol\beta_0'\mathbf{X})^2\}/\mathbb{E}\{(\boldsymbol\gamma_0'\mathbf{X})^4\}} = 3$, a reasonable balance between linear and quadratic signal parts. Despite the true $m_0(\cdot)$ being non-linear, we employed linear working models to update the nuisance posterior $\Pi_{\mathbf{m}}$ like {\tt Bols, Bridge} and {\tt Bsparse} methods as detailed in Section (ref). Note that due to potential misspecification, $\Pi_{\mathbf{m}}$ now contracts around a non-random limiting function $m^*(\mathbf{X}) := \widetilde\mathbf{X}'\boldsymbol\beta^*$, where $\widetilde\mathbf{X} = (1, \mathbf{X}')'$ and $\boldsymbol\beta^* := \arg \min_{\boldsymbol\beta \in \mathbb{R}^{p+1}} \mathbb{E} ( Y - \widetilde\mathbf{X}' \boldsymbol\beta )^2$, i.e., $\boldsymbol\beta^* = \{\mathbb{E}(\widetilde\mathbf{X}\widetilde\mathbf{X}')\}^{-1}\mathbb{E}(\widetilde\mathbf{X} Y)$, refer to Remark (ref).

Unlike the settings in Section (ref), where the theoretical ORE is attainable, it is {\it not} achievable here due to model misspecification. Instead, we calculated the {\it achievable} oracle asymptotic RE $(\mathrm{ORE}^*)$, defined as $\mathrm{ORE}^* := \mathrm{Var}(\widehat\theta_{\sup})/\tau^2_{n, N}(m^*)$, where $\tau^2_{n, N}(m^*) = \mathrm{Var}\{Y - m^*(\mathbf{X})\}/n + \mathrm{Var}\{m^*(\mathbf{X})\}/N$ and $m^*(\cdot)$ is the possibly misspecified limit of $\Pi_{\mathbf{m}}$. However, both $\mathrm{ORE}$ and $\mathrm{ORE}^*$ are reported as performance benchmarks.

Table (ref) and Tables (ref)--(ref) present the results on estimation and inference, respectively. Table (ref) shows that the REs of the SS estimators $\widehat\theta_{\mathrm{BDM}}$ and $\widehat\theta_{\mathrm{hBDM}}$, compared to the supervised estimator $\widehat\theta_{\sup}$, are substantially greater than 1, ranging roughly from 2.4 to 2.8 (matching the $\mathrm{ORE}^*$ closely) across most settings. Further, Tables (ref)--(ref) show that BDMI consistently achieves CovPs close to the nominal 95% level, and with significantly tighter CIs (typically 25--40% tighter) compared to the supervised approach across all settings and methods for $\Pi_{\mathbf{m}}$. All these findings highlight: (i) the {\it efficiency improvement} and (ii) {\it global robustness} that BDMI {continues} to enjoy {\it even under misspecification of $\Pi_{\mathbf{m}}$}, and further reinforces its first-order {\it insensitivity} to nuisance estimation bias. For more visual illustrations, see Figures (ref)--(ref) in the \hyperref[sec:supplementary]{Supplementary Material}.

A notable aspect of the RE values in Table (ref) is that the extent of the efficiency improvement is fairly {\it uniform} across the settings, and quite close to the {\it achievable} $\mathrm{ORE}^*$ in most cases -- with a slight lowering, in general, for the higher $p = 50$ case, as expected. This indicates no substantial additional finite sample losses in estimating $m^*(\cdot)$ under the low/moderate dimensional settings here. On the other hand, the difference between the achievable $\mathrm{ORE}^*$ and the optimal $\mathrm{ORE}$ indicate the (unrecoverable) difference due to the $O(1)$ bias stemming from $\Pi_{\mathbf{m}}$ targeting $m^*(\cdot)$ and not the true $m_0(\cdot)$. One notable exception to the general uniform pattern in the REs is the case of Bsparse for $p = s = 50$, where the REs, while still high, are closer to 2. This arises since the dense setting introduces an additional layer of (soft) misspecification, making consistent estimation of even the $m^*(\cdot)$ more challenging for a sparsity-friendly method at such a choice of $(p,s,n)$. Conversely, Bols and Bridge, which do not depend on sparsity, continue to provide higher REs around 2.5.

Finally, consistent with our findings in Section (ref), BDMI and h-BDMI {\it still} perform similarly across all settings, with h-BDMI having slightly higher REs, while also exhibiting some conservativeness in CovPs, at least in some cases. Furthermore, similar to Section (ref), the results (both for estimation and inference) remain fairly robust across $K = 5$ and $K=10$. Thus, we continue to recommend either choice in practice.

Overall, as shown in Sections (ref)--(ref), BDMI {\it always} achieves significant efficiency improvements and valid inference, under both correctly specified and misspecified models, thus validating our theoretical results.

Real data analysis: Application to NHEFS data

In this section, we apply the proposed BDMI approach to a subset of data from the National Health and Nutrition Examination Survey Epidemiologic Follow-up Study (NHEFS), a longitudinal study jointly initiated by the National Center for Health Statistics and the National Institute on Aging in collaboration with other agencies of the United States Public Health Service Hernan2020. The NHEFS was designed to investigate the effects of clinical, nutritional, demographic, and behavioral factors on various health outcomes, including morbidity and mortality. Data were collected during a baseline visit in 1971 and a follow-up visit in 1982. For our analysis, we focus on a cohort of 1425 individuals from this study. A detailed description of the dataset is available at \hyperlink{https://hsph.harvard.edu/miguel-hernan/causal-inference-book}{https://hsph.harvard.edu/miguel-hernan/causal-inference-book}. This dataset has been widely used in other studies for different purposes. For instance, Ertefaie2022 used this dataset to estimate causal parameters like the average treatment effect of quitting smoking on weight gain, and chakrabortty2022semi focused on quantile estimation under an SS framework.

Our primary goal is to estimate the {\it mean body weight}, $\theta_0$, of the entire cohort in 1982 under a semi-supervised framework. Additionally, we aim to investigate whether there is a significant change in body weight within the cohort between 1971 and 1982. To achieve this, we compared the analysis results to the baseline measurement from 1971, which had a mean of 70.99 and a standard error of 0.41 for the 1425 individuals. To benchmark our results, we also consider a {\it gold standard} scenario where all 1425 observations (for the response) are available in 1982 (which is the case for this data). We take the mean weight $\widehat \theta_{\mathrm{GS}} = 73.6$ of all 1425 individuals in the 1982 cohort as the gold standard ($\mathrm{GS}$) estimator (i.e., a close `proxy' of the truth).

To evaluate the performance of the proposed BDMI approach, we randomly select $n = 200$ observations as the labeled dataset $\mathcal{L}$, where body weight (response variable) is observed. For the unlabeled data $\mathcal{U}$, we randomly designate $N \in \{ 400, 800, 1220\}$ observations from the remaining data. This setup allows us to explore and compare the performance of BDMI under varying ratios of labeled and unlabeled data ($n/N$), particularly as this ratio approaches 0. In addition to body weight as the response variable, we considered a set of 20 important covariates in our analysis, including demographic, clinical, and behavioral factors (see Table (ref) in the \hyperref[sec:supplementary]{Supplementary Material} for their names and descriptions). These variables were also considered in other studies on this dataset, e.g., chakrabortty2022semi used them for SS quantile estimation.

The gold standard estimator $\widehat \theta_{\mathrm{GS}}$ provides a benchmark for evaluating and comparing the performance of BDMI versus the supervised approach (based on the labeled data only). Given the labeled and unlabeled data, we calculated the supervised posteriors $\Pi_{\sup}$ (see Section (ref)) and $\Pi_{\boldsymbol \theta} \equiv \Pi_{\mathrm{BDM}}$ based on the BDMI-CF approach (Algorithm (ref)) with $K =5$. From each posterior distribution, 1000 samples were obtained to compute the point estimators $\widehat \theta_{\sup}$ and $\widehat \theta_{\mathrm{BDM}}$, along with the respective $95\%$ credible intervals (CIs). Additionally, we calculated the {\it ratio} of the lengths ({\bf RL}) of the $95\%$ CIs from the supervised approach to those from BDMI. This RL serves as a natural measure of the relative efficiency of BDMI, where an RL greater than 1 indicates that BDMI provides tighter (and hence more efficient) CIs. Similar to our simulation study, 3 different methods are used to update the posterior $\Pi_{\mathbf{m}}$ of the nuisance parameter $m$, resulting in 3 distinct posteriors $\Pi_{\boldsymbol \theta} \equiv \Pi_{\mathrm{BDM}}$ for the BDMI approach. Table (ref) summarizes our findings from the data analysis.

Table (ref) highlights that BDMI demonstrates two key advantages over the supervised approach: {\it improved} accuracy and efficiency. First, the SS point estimates based on BDMI (across all versions) are consistently closer to the gold standard estimate, $\widehat \theta_{\mathrm{GS}} = 73.6$, compared to the supervised estimate $\widehat \theta_{\sup} = 72.6$ for all settings of $N$. Second, BDMI (across all versions) consistently produces significantly tighter $95\%$ CIs than the supervised approach, with efficiency gains quantified by the ratio of CI lengths (RL), ranging from $1.2$ to $1.7$ across all settings. This corresponds to $20-70\%$ tighter intervals for BDMI. Notably, for a fixed number of labeled data, as the ratio $n/N$ decreases (i.e., increasing the size of the unlabeled data), BDMI achieves substantial efficiency improvements by further reducing CI lengths compared to the supervised approach. For example, with $n=200$, increasing $N$ from 400 to 1220 improves the RL from around 1.24 to 1.74, reflecting a $40\%$ reduction in CI length for BDMI. These results indicate that the posterior spread under BDMI becomes increasingly tighter as more unlabeled data are incorporated. Hence, these findings highlight BDMI's ability to {\it efficiently} leverage unlabeled data, providing strong empirical support for our theoretical framework regarding the importance of $\lim_{n, N \to \infty} n/N = c \in [0,1)$. These results show that the BDMI procedure delivers both accurate point estimates (near identical to the GS version) and enhanced efficiency through shorter/tighter credible intervals, underscoring its advantage over the supervised approach. Finally, a notable feature of the BDMI based CIs for the mean weight of the 1982 cohort is that they consistently {\it exclude} the 1971 mean weight (70.99), indicating a significant {weight gain}, likely due to aging or quitting smoking Ertefaie2022. In contrast, the supervised approach {\it fails to detect} this change, as its $95\%$ CI $(70.6, 75.0)$ includes the 1971 mean. These results highlight the improved efficiency and higher {\it power} of BDMI for detecting significant (and scientifically meaningful) differences in weights between the two cohorts.

Concluding discussions

We proposed the BDMI procedure for estimating the population mean $\theta_0 = \mathbb{E}(Y)$ under the SS setting. To the best of our knowledge, this is the first attempt to establish a Bayesian method that achieves desirable SS inference goals, including {\it efficiency improvement} and {\it global robustness}, while providing {\it rigorous} theoretical guarantees. Our methodology ensures that the posterior $\Pi_{\boldsymbol \theta}$ of the parameter of interest $\theta$ contracts around the true parameter $\theta_0$ at the parametric rate $n^{-1/2}$ and is asymptotically Normal, {\it regardless} of the choice of method used to obtain a posterior for the nuisance parameter $m$, its contraction rate, or even potential misspecification of $m$. Moreover, the posterior mean of $\Pi_{\boldsymbol \theta}$, as an SS estimator of $\theta_0$, {\it always} possesses $\sqrt{n}$-consistency, asymptotic normality, and first-order {\it insensitivity}, in addition to being at least as {\it efficient} as the supervised estimator (sample mean of $Y$). These theoretical results have been rigorously established in Section (ref). One of the key contributions of BDMI lies in its ability to disentangle nuisance parameter estimation from inference on the target parameter by developing a novel debiasing approach under the Bayesian paradigm. It enables joint learning of the nuisance bias and the main parameter through targeted modeling of summary statistics, along with careful usage of sample splitting. We hope this research brings attention to the rarely used idea of modeling summary statistics within Bayesian inference and demonstrates its potential to address other Bayesian semi-parametric inference problems. While this work focuses on SS mean estimation, the underlying principles of BDMI can be extended to a broad range of problems, including missing data analysis, causal inference, and SS inference for other functionals. For instance, BDMI could be adapted to handle selection bias or distribution shifts between labeled and unlabeled data; this was recently explored in the frequentist SS literature zhang2021double but not yet addressed within a Bayesian framework. Further, extending BDMI to Bayesian SS inference for high dimensional target parameters (e.g., regression coefficients) poses additional theoretical and computational challenges, but also represents an important direction for future research. Finally, adapting BDMI’s debiasing framework to causal inference or missing data settings offers exciting opportunities for advancing Bayesian semi-parametric methodologies. We hope this work generates interest in considering such Bayesian problems in the future.