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.
81,096 characters · 10 sections · 65 citation commands
Bayesian Double Machine Learning for Causal Inference
\thispagestyle{empty}
\onehalfspacing
Two decades into the “credibility revolution” popularized by Mostly Harmless Econometrics, much work in applied microeconomics continues to rely on relatively simple reduced-form models to untangle cause-and-effect from observational data. In the simplest case, and the one we consider below, a researcher argues that a treatment $D$ is as good as randomly assigned after adjusting for a set of control variables $X$ and goes on to report the coefficient on $D$ from a regression of $Y$ on $(D,X)$ as the causal effect of interest. Naturally this selection-on-observables approach is only as good as the choice of controls; researchers who rely on it typically devote much of their paper to defending the choice of $X$ in their application.
All else equal, a selection-on-observables strategy is more convincing when it includes more observables. For this reason the dimension of $X$ can be very large relative to sample size. Even if there are relatively few “underlying” control variables, fully adjusting for them may require us to include a large number of interactions and other transformations in our regression. Applied researchers thus face a difficult trade-off: including more controls reduces the risk of bias but can lead to extremely imprecise estimates of the causal effect of interest.
In this paper, we provide a simple, fully Bayesian approach to navigating this bias-variance trade-off in a workhorse partially linear model: Bayesian Double Machine Learning (BDML). While our methods are Bayesian, the resulting point and interval estimators are designed to appeal to a broad audience. For Bayesians, BDML offers a natural way to bring the full Bayesian toolbox to bear on high-dimensional causal inference problems in applied microeconomics. For Frequentists, BDML provides superior finite-sample performance in practice--lower root mean-squared error and improved coverage of confidence intervals--while delivering the same attractive asymptotic guarantees as existing Frequentist approaches.
Any attempt to tame the estimation uncertainty that arises from including many controls relies on regularization, also known as shrinkage. Rather than “projecting out” $X$ from $Y$ and $D$ using ordinary least squares (OLS), a regularized estimator constrains the high-dimensional vector $\beta$ of coefficients on $X$ in the regression of $Y$ on $(D,X)$. Bayesian approaches do so by placing a prior on $\beta$; machine learning (ML) approaches do so either by adding a penalty term to the objective function to favor “smaller” values for the elements of $\beta$, as in LASSO or Ridge, or by penalizing model complexity through implicit algorithmic constraints, as in random forests and gradient boosting.
Whether Bayesian or Frequentist, implicit or explicit, any application of regularization to causal inference problems faces a key challenge: regularization-induced confounding (RIC).\footnote{This term was coined by hahn_regularization_2018.} Unlike OLS, regularized estimators only partially adjust $D$ and $Y$ for the observed confounders $X$. This induces a correlation between $D$ and the regularized regression residuals, biasing the estimated causal effect. Granted, researchers must be prepared to accept at least some increase in bias in exchange for reduced estimation uncertainty. The problem of RIC is that this bias can be extremely sensitive to the precise way in which regularization is carried out, depending on the relationship between $D$ and $X$. Indeed, na\"{i}ve approaches to regularized estimation can easily increase the root mean-squared error (RMSE) of our causal effect estimator, adding a bias that far outweighs the reduction in variance. Crucially, our BDML approach side-steps this problem completely: using a simple re-parameterization of the model, we avoid RIC while relying on standard methods for Bayesian linear regression.
From a subjective Bayesian point of view, RIC may seem puzzling. If we believe our prior and likelihood, Bayes optimality requires us to base estimation and inference on the implied posterior. In practice, however, priors are often chosen for pragmatic reasons, and as sims_robins-wasserman_2012 notes, such priors “can unintentionally imply dogmatic beliefs about parameters of interest” when expanded “unthinkingly to high dimensions.” linero_nonparametric_2023 shows that the causal inference problem we study here has this problem in spades. It is common, and convenient in practice, to place independent priors on $\gamma$ and $\beta$, the coefficients on $X$ in the “propensity score” and “outcome” regressions, respectively. In low dimensions, this is innocuous. In high dimensions, on the other hand, it implies a dogmatic prior belief that there is no selection bias. If we truly believed this, of course, there would be no reason to adjust for $X$ in the first place!
To avoid prior dogmatism without the need to elicit a dependent prior for $(\gamma, \beta)$, BDML replaces the outcome regression of $Y$ on $(D, X)$ with a reduced form regression of $Y$ on $X$. This leaves us with a bivariate, seemingly unrelated regression of $(Y,D)$ on $X$. Under the selection-on-observables assumption, the error terms of these two equations are correlated, and their covariance matrix contains all the information needed to carry out estimation and inference for the causal effect of $D$ on $Y$. We show that placing independent priors on the coefficients of the seemingly unrelated regression equations runs no risk of prior dogmatism, regardless of the dimension of $X$. Any standard prior on the error covariance matrix implies a plausible range for the a priori selection bias.
We call our method BDML to highlight its connection to the Frequentist Double Machine Learning (FDML) approach that has become extremely popular in recent years chernozhukov_doubledebiased_2018. Like our proposal, FDML can be viewed as a means of avoiding RIC when carrying out causal inference in high dimensions. In the partially linear regression model we study here, FDML likewise considers a bivariate reduced form regression of $(Y,D)$ on $X$. After estimating the coefficients of these two regressions using an ML approach, FDML recovers the causal effect of interest using a further regression of the $Y$-residuals on the $D$-residuals. This approach is computationally convenient and can be given a rigorous large-sample justification. In practice, however, it amounts to treating the reduced form coefficients as if they were known--that is, FDML profiles out the nuisance parameters by plugging in point estimates. When the dimension of $X$ is large relative to sample size there may be considerable sampling uncertainty in the reduced form coefficients. BDML explicitly acknowledges this by marginalizing over all sources of uncertainty in the problem, leading to improved small-sample performance as measured by traditional Frequentist metrics, e.g.\ root mean-squared error (RMSE) and coverage probability.
Crucially, the finite-sample advantages of BDML do not come at the expense of large sample guarantees. To show this, we study a fixed coefficient asymptotic sequence in which the number of control regressors $p$ grows with sample size $n$, matching the approach of the FDML literature. We consider three estimators: BDML, FDML, and a “n\"{a}ive” Bayesian alterantive that ignores RIC. Under regularity conditions, we show that all three yield consistent estimators of the causal effect of interest. When they are $\sqrt{n}$-consistency, we likewise show that all three have the same asymptotic variance. The asymptotic bias of BDML and FDML, however, is strictly smaller than that of the n\"{a}ive Bayesian approach: $p^2/n^2$ compared to $p/n$. Moreover, the n\"{a}ive estimator is only $\sqrt{n}$-consistent if the number of controls is much smaller than sample size: $p/\sqrt{n} \rightarrow 0$. BDML and FDML, in contrast, allow for a much larger number of controls, requiring only $p/(n^{3/4})\rightarrow 0$. Asymptotically, BDML improves on the n\"{a}ive Bayesian approach while matching the performance of FDML.
Continuing our asymptotic analysis, we build on results from walker_parametrization_2025 to establish a Bernstein-von Mises (BvM) theorem for BDML when $(\log n)^2 \prec p \prec \sqrt{n}$.\footnote{We use the symbol $a \prec b$ to indicate that $\lim_{n\rightarrow \infty} a/b \rightarrow 0$.} This result further strengthens the Frequentist appeal of BDML, showing that posterior credible sets constructed using our approach are valid Frequentist confidence intervals. Our BvM result further establishes the semi-parametric efficiency of BDML while imposing weaker assumptions on the prior than the recent Bayesian literature luo_semiparametric_2023,breunig_semiparametric_2024. While we prove the BvM theorem under the assumption of Gaussian errors, these results are robust to mis-specification of the error distribution.
We conclude by comparing the finite-sample performance of BDML against a range of Bayesian and Frequentist competitors in a comprehensive simulation study and a number of empirical examples drawn from the recent literature in applied microeconomics. Our results show that BDML provides the “best of both worlds”, matching the asymptotic properties of Frequentist alternatives while delivering dramatically improved finite-sample performance. Compared existing Bayesian alternatives hahn_bayesian_2020,linero_nonparametric_2023, BDML maintains its impressive performance across a much wider range of data-generating processes (DGPs). Since BDML is fully-Bayesian, it correctly accounts for all sources of sampling uncertainty, smoothly adapting to key features of the DGP.
This paper relates to a large recent literature on (Frequentist) double machine learning (FDML), also known as “double-debiased machine learning”. See ahrens_introduction_2025 for a comprehensive practical introduction and chernozhukov_doubledebiased_2018 for further theoretical details. FDML provides a general framework for estimating a low-dimensional target parameter using moment conditions that depend on a high-dimensional nuisance parameter. The na\"ive plug-in approach--constructing a preliminary ML estimate of the nuisance parameter and plugging it into the original moment condition--has two serious shortcomings. The first is regularization bias, a more general term for RIC. The second is overfitting bias from using the same data twice--first to estimate the nuisance parameter and then to estimate the target parameter. By orthogonalizing the second-step moment condition with respect to the nuisance parameter, FDML inoculates against regularization bias to first order. By using sample splitting, it likewise avoids overfitting bias.
Two key advantages of FDML are its simplicity and generality: whenever we can derive the appropriate orthogonalized moment condition, FDML justifies treating the estimated nuisance parameter as if it were known while providing theoretical guarantees across a wide range of first-step ML estimators. These attractive properties have led applied researchers in economics to employ FDML in a range of empirical setting, especially the “post-double selection” (PDS) LASSO approach of belloni_inference_2014. Recent papers that apply FDML as their primary method in an observational setting include appl_dube_monopsony_2020 who estimate monopsony power in online labor markets, and appl_yang_bign_audit_2020 who estimate the “Big N” audit quality effect. While its primary motivation is to reduce bias in observational settings, FDML is also widely applied in experimental settings to flexibly incorporate covariates, either as a robustness check appl_alsan_representation_clintrials_2023, appl_armand_resource_curse_2020 or with the aim of improving estimation precision appl_hussam_psych_refugees_2022,appl_cotti_ecig_flavor_bans_2025,appl_delfino_pink_collar_2024, appl_chioda_skillstraining_uganda_2021. FDML has also been used to evaluate, and ultimately question, the plausibility of the selection-on-observables assumption, either by comparing against an instrumental variables (IV) benchmark appl_duflo_ghana_2017 or by creating a “simulated” observational dataset from a collection of randomized controlled trials appl_gordon_ms_advertising_measure_2023.
Although FDML is indisputably a valuable tool in applied economics, some recent work suggests that its finite sample performance can be sensitive to details of the first-step ML procedure in ways that may lead it to underperform in practice angrist_machine_2022, wuthrich_omitted_2023, bach_hyperparameter_2024. ahrens_introduction_2025 in particular stress that “poorly tuned or ill-suited nuisance estimators can produce severely misleading point estimates”. As such, they encourage researchers using FDML to “carefully choose and tune their nuisance function estimators and consider a rich set of estimation approaches.” While by no means an insurmountable challenge, the need for careful tuning makes FDML more complicated to implement in practice than its theoretical simplicity may suggest. In our view, this challenge makes Bayesian alternatives to FDML particularly appealing. By marginalizing over uncertainty in the nuisance parameter, BDML provides a fuller accounting for finite-sample estimation uncertainty. Moreover, our method can provide robust performance across a wider range of data generating processes through the use of hierarchical priors, which adapt the degree of shrinkage to the data at hand.
Our work also relates to a growing literature on Bayesian approaches to causal inference: see heckman_treatment_2014, li_bayesian_2023, and linero_how_2023 for broad discussions of this field. The studies closest to ours are hahn_regularization_2018,hahn_bayesian_2020 and linero_nonparametric_2023, who consider the same partially linear model that we analyze below.\footnote{See belloni_inference_2014 for a frequentist analysis of the partially linear model.} Compared to BDML, which marginalizes over all sources of uncertainty, these papers adopt a “hybrid” strategy: profiling over some parts of the problem and marginalizing over others. This hybrid approach is computationally convenient and can be given a sound Bayesian justification. Our simulation studies reveal, however, that the full marginalization provided by BDML is crucial for obtaining inferences that have good frequentist properties across a wide range of data generating properties. As mentioned above, our BvM results draw heavily on walker_parametrization_2025. Recent Bayesian work exploring issues related to those considered here, albeit in different models, includes saarela_bayesian_2016, ray_semiparametric_2020, antonelli_causal_2022,luo_semiparametric_2023, breunig_semiparametric_2024 and breunig_double_2025.
The remainder of this paper is organized as follows. (ref) introduces our model and notation while (ref) discusses the problem of RIC in detail. (ref) presents our Bayesian Double Machine Learning approach and describes how it relates to both Bayesian and Frequentist alternatives, while (ref) presents a comparison of the asymptotic properties of BDML, FDML, and a na\"{i}ve Bayesian alternative. Finally, (ref) presents a simulation study comparing the finite-sample performance of BDML against a range of Bayesian and Frequentist competitors. Proofs appear in (ref).
We consider the problem of estimating and carrying out inference for the parameter $\alpha$ in a partially linear structural model of the form
where $Y_i$ is the outcome of interest, $D_i$ is the “treatment” or “policy” variable, and $X_i$ is a vector of control variables.\footnote{Our model and parameter of interest match those of belloni_inference_2014.} We do not assume that $D_i$ is binary; if it is, (ref) can be viewed as a special case of the “selection on observables” framework. For reasons that will become clear below, we introduce an additional reduced form model, namely
where $m(X_i) \equiv \mathbbm{E}[D_i|X_i]$ so that (ref) holds by construction and we have
by iterated expectations and $\mathbbm{E}[\varepsilon|X_i,D_i]=0$, since $\varepsilon_i$ and $V_i$ are mean zero.
A leading example of (ref)--(ref) assumes that $g(X_i) = X_i'\beta$ and $m(X_i)=X_i'\gamma$ where $\beta$ and $\gamma$ are $p$-vectors of unknown regression coefficients. This simplification gives
where we have
While the methods we develop in this paper are easily extended to the general case of (ref)--(ref), we focus on the special case of (ref)--(ref) to simplify the exposition. We assume throughout that the dimension $p$ of $X_i$ is large relative to sample size $n$. If the dimension of $X_i$ exceeds the number of observations, then there is no unique ordinary least squares (OLS) estimator of $\alpha$ based on (ref). Even if there are more observations than regressors, the OLS estimator of $\alpha$ is likely to be extremely noisy if $X_i$ if contains many “relatively unimportant” control variables--variables that, taken individually, have little predictive power for $D_i$ and $Y_i$. This situation is arguably quite common in practice. Applied researchers may know which controls it would suffice to adjust for without being sure of which are most important.
As is well-known, the OLS estimator for $\alpha$ corresponds to the Bayesian posterior mean under (ref) with a normal distribution for $\varepsilon_i$ and flat priors on $(\alpha, \beta)$. Accordingly, in a situation where we expect OLS to perform poorly because the dimension of $X_i$ is large, a natural idea would be to “shrink” our estimates by introducing an informative prior on $\beta$. While there are many possibilities for the choice of shrinkage prior, the message of this section is that none of them is likely to perform well in practice. To illustrate this point, suppose we retain a flat prior for $\alpha$ but place independent normal priors on the elements of $\beta$:
Under the simplifying assumption that $\sigma^2_\varepsilon$ is known, the Bayes estimator for $(\alpha, \beta)$ under squared error loss has a closed-form solution:
As seen from (ref), the Bayes estimator in this setting corresponds to a modified version of Ridge Regression in which $\alpha$ is unpenalized. As usual in Ridge Regression, we implicitly assume that $D_i$, $X_i$ and $Y_i$ have all been de-meaned so that the regression requires no intercept. If (ref) represents genuine subjective researcher beliefs over $\beta$ then, from a Bayesian point of view, $\widehat{\alpha}_\lambda$ is the mean-squared error optimal estimator of $\alpha$. In practice, however, shrinkage priors like (ref) are often used for pragmatic reasons: when $X_i$ is high-dimensional, it may be difficult if not impossible to elicit a fully-informative prior over $\beta$. For the purposes of predicting $Y_i$ this is usually innocuous, provided that $\sigma_\beta^2$ is chosen in a reasonable way. But when the goal is to understand the causal effect of $D_i$, the situation changes dramatically. To illustrate this point, we derive the (Frequentist) bias and variance of $\widehat{\alpha}_\lambda$ as follows under the simplifying assumption that $\varepsilon$ is homoskedastic.\footnote{hahn_regularization_2018 derive a special case of $\text{Bias}(\widehat{\alpha}_\lambda)$, but their expression contains a minor error.}
The $j$\textsuperscript{th} element of $\widehat{\rho}$ from (ref) is the slope coefficient from a regression of $X_i^{(j)}$ on $D_i$ without an intercept (both have been de-meaned) while the $j$\textsuperscript{th} column of $\widehat{\xi}$ is the corresponding vector of residuals. Setting $\lambda = 0$ corresponds to $\sigma_\beta^2 \rightarrow \infty$ and gives the OLS estimator of $\alpha$. In this case the bias is zero and the variance simplifies to \[ \text{Var}(\widehat{\alpha}_{\text{OLS}}) = \sigma^2_\varepsilon \left[ (D'D)^{-1} + \widehat{\rho}'R^{-1} \widehat{\rho}\right]. \] If $\sigma^2_\beta$ is finite, then $\lambda \neq 0$ and $\widehat{\alpha}_\lambda$ is biased unless $\widehat{\rho}$ or $\beta$ equals zero. Of course this in and of itself is not a convincing argument against shrinkage. In general, higher values of $\lambda$ imply lower values of $\text{Var}(\widehat{\alpha})$. Even if $n > p$ so that OLS remains feasible, we should still prefer to set $\lambda > 0$ provided that the decrease in variance outweighs the increase in (squared) bias.
The problem isn't the presence of bias but rather its dependence on other aspects of the model. We include covariates on the right-hand side of (ref) precisely because we believe that $\alpha$ would lack a causal interpretation without them. This implies an a priori belief that the control variables covary with both the outcome and treatment. As seen from (ref), for a fixed value of $\lambda$ the bias of $\widehat{\alpha}_\lambda$ involves both $\beta$ and $\widehat{\rho}$. All else equal, larger values of $\widehat{\rho}^{(j)}\beta^{(j)}$ imply a higher bias. Our beliefs about the magnitude of $\beta^{(j)}$ are already accounted for, at least in principle. Since $\lambda = \sigma_\varepsilon^2 / \sigma_\beta^2$, a belief that $|\beta^{(j)}|$ will be large on average leads us to shrink less, using a smaller value of $\lambda$ and keeping the bias in check. But $\text{Bias}(\widehat{\alpha}_\lambda)$ also depends crucially on $\widehat{\rho}$, and $\lambda$ fails to take this into account.\footnote{The stated priors on $(\alpha, \beta)$ along with the fixed known value of $\sigma_\varepsilon^2$ implicitly encode beliefs about the likely magnitude of $\widehat{\rho}$, as we explain further below.} If the treatment and controls are strongly correlated, even a reasonable choice for $\sigma_\beta^2$ can produce a highly biased estimator of $\alpha$, a phenomenon that hahn_regularization_2018 call regularization-induced confounding (RIC). The name alludes to the fact that “regularizing” the OLS estimator by setting $\lambda > 0$ causes the residuals to violate the sample analogue of $\mathbbm{E}[D_i\varepsilon_i|X_i]=0$, the population moment condition that gives $\alpha$ a causal interpretation. While RIC is typically equated with unacceptably high bias, numerical experiments using (ref)--(ref) confirm that RIC implies an unfavorable bias-variance trade-off. This is the sense in which we use the term.
The trouble, it seems, is failing to account for our prior beliefs about $\gamma$, the reduced form coefficient from (ref) that reflects the dependence between $D_i$ and $X_i$. A natural solution, then, would be to incorporate these beliefs into our model. Doing so requires us to add an equation and work with a bivariate regression including both (ref) and (ref). Staying as close as possible to the single-equation example from above, suppose we specify a normal likelihood for $(\varepsilon_i, V_i)$ and place independent normal shrinkage priors on the elements of $\beta$ and $\gamma$, all centered at zero. In this model, an additional hyperparameter $\sigma_\gamma^2$ quantifies our prior beliefs about the strength of dependence between the treatment and control variables. Unfortunately, $\sigma_\gamma^2$ will have no effect whatsoever on our estimates or inferences for $\alpha$. Recall from (ref) above that the causal assumptions of our model require $\text{Cov}(\varepsilon_i, V_i) = 0$. Adopting the shorthand $\theta$ to refer to all model parameters, this implies that the likelihood factorizes as follows
Notice that the first factor involves only $(\alpha, \beta, \sigma_\varepsilon^2)$ while the second involves only $(\gamma, \sigma_V^2)$. Therefore, unless we specify a prior in which $\beta$ and $\gamma$ Treating $\sigma_\varepsilon^2$ as known, the bias and variance expressions from (ref)--(ref) continue to hold and the problem of RIC remains. linero_nonparametric_2023 calls a causal model in which the likelihood factorizes as in (ref) and $\beta$ and $\gamma$ are a priori independent a Bayesian-ignorable model, and points out that such models imply a highly dogmatic prior for the extent of selection bias when $p$ is large. This is another lens through which to view RIC: it arises from “accidentally” imposing a strong prior belief that there is little or no selection bias. We explore this phenomenon in more detail in (ref) below.
In principle, specifying an informative prior that incorporates dependence between $\beta$ from (ref) and $\gamma$ from (ref) would solve the problem of RIC introduced in (ref). Eliciting such a prior in practice, however, may be extremely challenging, especially when the number of control variables is large. This motivates us to introduce an alternative approach, one that is inspired by ideas from the recent Frequentist literature on so-called “double-debiased machine learning”, often referred to as “double machine learning” for short. Substituting (ref) into (ref) gives the reduced form regression
where $U_i \equiv \varepsilon_i + \alpha V_i$ and $\delta \equiv \alpha \gamma + \beta$. Since since $\mathbbm{E}[\varepsilon_i|X_i,D_i]=0$ by (ref) and $X_i$ is uncorrelated with $V_i$ by (ref), $X_i$ is uncorrelated with $U_i$. Assuming normal errors for simplicity, this gives the following bivariate reduced form regression model
Given that our goal is to estimate and carry out inference for $\alpha$, it may seem odd that we have eliminated this parameter from the regression equation for $Y$. Crucially, however, the errors $(U_i, V_i)$ in (ref) are correlated whenever $\alpha\neq 0$. This is what will allow us to learn the causal effect of interest. By (ref) we have $\text{Cov}(\varepsilon,V) = 0$ and thus
by the definition of $U$. Therefore,
Notice that knowledge of $\Sigma$ immediately implies knowledge of $\alpha$ via the simple transformation
The relationship in (ref) suggests the following approach to estimation and inference for $\alpha$.
The procedure given in (ref) is simple, flexible, and fully Bayesian. At the same time, it allows us to avoid the problem of inadvertently specifying a dogmatic prior on selection bias, pointed out by linero_nonparametric_2023 and explained in (ref) above. First, because we work with a pair of reduced form regressions rather than a structural regression for $Y_i$ and a reduced form regression for $D_i$, our likelihood does not factorize: $(U_i,V_i)$ are correlated whereas $(\varepsilon_i,V_i)$ are not. Second, because $\delta = (\alpha \gamma + \beta)$, placing independent priors on the reduced form regression coefficients $\delta$ and $\gamma$ does not imply independent priors for $\beta$ and $\gamma$. Hence, (ref) in general satisfies neither of the two conditions required for Bayesian-ignorability introduced by linero_nonparametric_2023. Naturally, the properties of point and interval estimates based on (ref) will depend on the specific choice of priors for $(\delta, \gamma, \Sigma)$. We explore this point in detail below and and compare the theoretical and practical performance of different priors. The key insight of our approach, however, is that re-writing the model from (ref)--(ref) in the form given by (ref) allows us to side-step the problem of RIC using a standard Bayesian multivariate regression model with conditionally conjugate priors. As far as we are aware, this observation is new to the literature.
Because the estimator $\widehat{\alpha}_\lambda$ from (ref) can be viewed from a Frequentist perspective as an instance of Ridge Regression, (ref)--(ref) imply that RIC is just as much of a challenge for frequentists who apply machine learning methods to high-dimensional causal inference problems as it is for Bayesians. The increasing recognition of this problem over the past decade has led to the development of a framework called (frequentist) double machine learning (FDML). See chernozhukov_doubledebiased_2018 for an overview. The simplest FDML approach to estimating $\alpha$ in (ref) proceeds by running three regressions. First, let $\hat{\delta}_{\text{ML}}$ be the estimated coefficient vector from some unspecified “machine learning” (i.e.\ regularized) regression of $Y_i$ on $X_i$, e.g.\ ridge regression, LASSO, random forests, etc.\ Second, let $\hat{\gamma}_{\text{ML}}$ be the estimated coefficient vector from some unspecified regularized regression of $D_i$ on $X_i$. The FDML estimator is then obtained by regressing the residuals from the $Y_i$ on $X_i$ regression on the residuals from the $D$ on $X$ regression, in other words
Like (ref) from above, (ref) involves reduced form regressions of $Y$ and $D$ on $X_i$. Because of this similarity, we refer to the approach from (ref) as Bayesian Double Machine Learning (BDML). From a Bayesian perspective, BDML provides a way to guard against RIC while using a fully generative model and conducting inference in a manner that respects the likelihood principle. This gives researchers access to the full arsenal of Bayesian tools, from hierarchical modeling to posterior predictive checks. But even for researchers who are primarily concerned with frequentist performance, BDML offers two key advantages. First, it allows practitioners to incorporate subject matter expertise into estimation through the choice of prior, potentially yielding sizeable efficiency gains in small samples. To show how this can be achieved in practice, we propose a simple and effective way of incorporating researcher beliefs concerning the R-squared of the reduced form regressions, building on the so-called “R2D2” prior of zhang_bayesian_2022. Second, rather than simply “plugging in” point estimates of $\delta$ and $\gamma$, as in FDML, our BDML approach marginalizes over these high-dimensional parameters, correctly accounting for all sources of estimation uncertainty for $\alpha$. As such, BDML may provide more accurate frequentist inference while avoiding the need for cumbersome and computationally-expensive “cross-fitting” procedures, as are commonly employed in applications of FDML.
While ours is not the only Bayesian proposal for avoiding RIC in a model like (ref), BDML differs in important ways from existing methods. The two closest to our approach are the methods of hahn_regularization_2018,hahn_bayesian_2020 on the one hand and linero_nonparametric_2023, on the other hand. Both begin with a preliminary Bayesian linear regression of $D_i$ on $X_i$ to yield a Bayes estimator of $\gamma$ that is then substituted into a second Bayesian regression.\footnote{While hahn_regularization_2018 takes a more involved joint estimation approach, their followup paper hahn_bayesian_2020 advocates the simpler two-step approach that we describe here.} The approaches differ in the precise regression model that they estimate in the second step. hahn_regularization_2018 rely on the same substitution as our (ref) but, rather than absorbing $\alpha V_i$ into a reduced form error term, replace it with its point estimate $\widehat{V}_i$ based on a first-step regression.
Note that the parameter $\delta_2=\delta+\alpha \widehat{\gamma} $ from (ref) does not in general coincide with our $\delta$ from (ref). linero_nonparametric_2023 generalizes the second-stage regression from hahn_regularization_2018. Defining $\widehat{D}_i \equiv X_i'\widehat{\gamma}$, we can re-express the regression from (ref) as \[ Y_i = \alpha D_i - \alpha \widehat{D}_i + X_i'\delta_2 + \nu_i. \] linero_nonparametric_2023 estimates the parameters of this regression without the restriction that the coefficient on $\widehat{D}_i$ equals the negative of the coefficient on $\widehat{D}_i$.
Notice that the regression from (ref) has an identification problem: $\widehat{D_i} = X_i\widehat{\gamma}$ is perfectly collinear with $X_i$. We write $\delta_3$ as the vector of regression coefficients for $X_i$ to emphasize that, when estimated without the restriction that $\kappa = -\alpha$, this procedure will not coincide with (ref).
Compared to the procedures in Algorithms (ref)--(ref), our BDML approach differs in two important ways. First, our inferences for $\alpha$ incorporate posterior uncertainty over both $\gamma$ and $\delta$. hahn_bayesian_2020 argue that, since Bayesian inference in a regression model conditions on the realizations of the regressors, a two-step approach can be viewed as fully-Bayesian: any Bayes-estimator of $\gamma$ from the first-step regression is necessarily a measurable function of $(D, X)$, so plugging this estimate into a second-step regression does not violate the likelihood principle. We agree with this argument, as far as it goes. But there are countless approaches that adhere to the likelihood principle, all of which yield different posteriors. As we show in our simulations below, accounting for uncertainty in $\widehat{\gamma}$ is crucial for obtaining inferences that have good frequentist properties in addition to a sound Bayesian justification. Second, by residualizing both $D_i$ and $Y_i$ with respect to $X_i$, BDML achieves robustness properties similar to those of frequentist double machine learning in a fully-Bayesian model. In contrast hahn_regularization_2018 residualizes only $D_i$. Alluding to similar robustness concerns raised by zigler_central_2016, linero_nonparametric_2023 suggests that his two-step approach could give added robustness when either (ref) or (ref) may be misspecified by eliminating “feedback” from the estimation problem. As we show below, however, our approach performs well in practice even in the simulation designs from linero_nonparametric_2023. This is because our approach, like FDML, is “double de-biased” where as Algorithms (ref)--(ref) are only “single de-biased”, a point we expand upon below.
We now explore the large-sample properties of the BDML approach from (ref) and contrast them with those of the frequentist alternative, FDML, as well as those of a “na\"{i}ve” Bayesian approach analogous to the Ridge Regression estimator from (ref). We first introduce some notation. Let $X$ be the $(n\times p)$ matrix of control regressors from above and define the $(n\times 2)$ matrix of “outcomes” $W$, the $(p\times 2)$ matrix of coefficients $B$, and the $(n \times 2)$ matrix of errors $E$ according to
where $E_i' \equiv (U_i, V_i)$ denotes a specified row of $E$. Then, we can express the reduced form regression model from (ref) in matrix form as
We consider an asymptotic sequence in which the number of control regressors $p$ grows with sample size $n$ and the true parameter values are fixed. We write these fixed true values as $\Sigma^*$ and $B^* \equiv
$ to distinguish them from the random variables $\Sigma$ and $B$ that represent prior (or posterior) uncertainty in a Bayesian model. Analogously, when we have occasion to refer to the true values of the structural parameters from \eqref{eq:linearY}, we write them as $(\alpha^*,\beta^*\sigma_\varepsilon^*)$.
Both for simplicity and to underscore the fact that BDML does not require the researcher to specify unusual or elaborate priors, we work with the classic conditionally conjugate prior for seemingly unrelated regression (SUR) models.\footnote{See zellner_introduction_1971 section 8.5 for a textbook treatment.} In particular, we combine (ref) with the prior $\pi(\Sigma,B) = \pi(\Sigma)\pi(\delta)\pi(\gamma)$ where
and $(\tau_\delta, \tau_\gamma)$ denote the prior precisions for the reduced-form parameters $\delta$ and $\gamma$.\footnote{As we explain further below, our use of a normal likelihood for $E_i$ is innocuous.} For purposes of comparison, we also consider a na\"ive Bayesian estimator given by the Ridge Regression estimator from (ref). This approximation treats the error variance $\sigma_\varepsilon^2$ as a known constant and imposes the prior $\beta \sim \text{Normal}_p(0, \mathbbm{I}_p/\tau_\beta)$, where $\lambda = \tau_\beta \sigma_\varepsilon^2$
Equation (ref) describes the Bayesian model we study. We now describe the assumptions that we place on the true data-generating process. For ease of notation, we state these in terms of the reduced form model. Combining these assumptions with (ref)--(ref) implies closely related assumptions for the structural model. The causal structure, in turn, gives a meaningful causal interpretation to $\alpha = \sigma_{UV}/\sigma_V^2$. We first assume that the researcher observes a sample of $n$ iid observations from the reduced form model.
While (ref) imposes normality of the errors $E_i$, as we explain further below this assumption turns out to be innocuous. Our next assumption concerns the true reduced form parameters $\Sigma^*$ and $(\delta^*, \gamma^*)$. As mentioned above, we consider an asymptotic sequence in which the number of control regressors, $p$, grows with sample size. This implies that the dimensions of $\delta$ and $\gamma$ likewise grow. To accommodate this, we consider $\delta^*$ and $\gamma^*$ to be square-summable infinite sequences of parameters, i.e.\ elements of the sequence space $\ell^2(\mathbbm{N})$. In finite samples only the first $p$ elements of each are “active”, but as we obtain more control regressors, the number of regression coefficients increases to match. Our assumption ensures that these coefficients do not “explode” as the number of regressors increases. Intuitively, each additional control matters “less on average” the more that we include in our model. Because the dimension of $\Sigma^*$ does not change with sample size, we need only assume that this variance-covariance matrix is well-behaved.
Our next assumption concerns the control regressors $X_i$. For simplicity, and without loss of generality, we assume that $X_i$ has been centered around its mean. To make it easier to compare our BDML procedure against its frequentist cousin, we further assume that there are more observations than controls: $p < n$. The remaining conditions we impose on $X_i$ restrict the variances and covariances of the control regressors as $p$ grows. To state them, let $\Sigma_X\equiv \text{Var}(X_i)$ and $\lambda_j(\Sigma_X)$ be the $j$-th eigenvalue of $\Sigma_X$, where these values are arranged in descending order. Further let $F_p(t) \equiv \frac{1}{p} \sum_{j=1}^{p} \mathbbm{1}\left\{\lambda_j(\Sigma_X) \leq t\right\}$ be the empirical spectral distribution of $\Sigma_X$.
Parts (iii) and (iv) of (ref) ensure that the spectrum of $\Sigma_X$ remains “well-behaved” as the dimension of $X_i$ grows. The final ingredient of our asymptotic framework is an assumption about the relative rates of $n$, $p$, and the prior precisions. To state it we use the notation $a_n \asymp b_n$ to indicate that two sequences are “of the same order”, i.e.\ that $a_n = O(b_n)$ and $b_n = O(a_n)$.
Part (i) of (ref) states that the sample size dominates the number of controls and is the asymptotic equivalent of $p < n$ from (ref) (i) above.\footnote{Formally (ref) (i) is redundant given (ii) and (iii), but we state it explicitly to emphasize a key feature of our asymptotic framework.} Since $p \rightarrow \infty$, this means that there are many controls but not so many that they overwhelm the available data. Part (ii) of (ref) states that the sample size dominates each of the prior precisions. Intuitively, this means that the priors are only “weakly informative” in the sense that they will be dominated by the data in large samples. Along this asymptotic sequence, the shrinkage bias of the Ridge estimator from (ref) vanishes in the limit. An $o(n)$ rate for the prior precisions also accords with the typical assumption from the frequentist literature on double machine learning. Finally, part (iii) asserts that the prior precisions are of the same order as the number of controls. This means that our prior imposes more shrinkage when there are more controls. Taken together, assumptions (ref)--(ref) implied that, under the prior from (ref), the implied R-squared values from the regressions of $D$ on $X$ and $Y$ on $X$ remain well-behaved in the limit: they do not diverge to zero or one. See (ref) in the Appendix for details.
Our BDML approach places priors on the reduced form coefficients $(\Sigma, B)$. We now consider the prior that (ref) implies for other quantities of interest. We begin with the causal effect of interest, $\alpha$. Partition $\Sigma$ into blocks $\Sigma_{11}, \Sigma_{12}$ and $\Sigma_{22}$ and likewise partition $\Sigma_0$ into $\Sigma_{0,11}$ and so on. Then $\alpha = \Sigma_{12}/\Sigma_{22}$. Under (ref), our prior for $\alpha$ follows a location-scale t-distribution.
If the prior scale matrix $\Sigma_0$ is diagonal, as is typical in practice, then $\Sigma_{0,12}/\Sigma_{0,22} = 0$ so the implied prior for $\alpha$ is symmetric around zero. Choosing a small value for $\nu_0$ gives this prior heavy tails, making it compatible with both very large positive and negative values for the effect of interest. This is reassuring: although we have placed a prior on the reduced form rather than the structural model, we can easily avoid imposing unintentionally strong prior beliefs on $\alpha$ by making $\Sigma_0$ diagonal.
We now consider the implied prior for a different but equally important quantity, namely the selection bias. This idea is most easily explained in the case where $D_i$ is binary. For the general version, see Propositions (ref)--(ref) in the Appendix. Following the usual convention in applied microeconomics, we define the selection bias as the difference between the quantity identified by a simple comparison of mean outcomes, treated minus untreated, and the true causal effect. Under (ref), this becomes \[ \text{SB} \equiv \left[ \mathbbm{E}(Y_i|D_i=1) - \mathbbm{E}(Y_i|D_i=0) \right] - \alpha = \left[ \mathbbm{E}(X_i|D_i=1) - \mathbbm{E}(X_i|D_i=0)' \right]\beta. \]
As mentioned in (ref) above, linero_nonparametric_2023 shows that any model that factorizes as in (ref) and in which and $\beta$ and $\gamma$ are a priori independent implies a highly dogmatic prior for the extent of selection bias when $p$ is large. This is another lens through which to view RIC: it arises from “accidentally” imposing a strong prior belief that there is little or no selection bias. The following result shows that the na\"ive estimator suffers from this problem whereas our BDML estimator does not.
When $n$ and $p$ are large, the na\"ive approach in effect assumes away selection bias: it implies a prior that is concentrated in a small neighborhood around zero. But if there is no selection bias, then there is no need to adjust for $X_i$ in the first place: it is as if $D_i$ had been randomly assigned! Intuitively, in high-dimensional spaces randomly chosen vectors tend to be nearly orthogonal to each other. Since $\Sigma_X$ is positive definite, the numerator \(\gamma'\Sigma_X\beta\) from $\text{SB}_{\text{naive}}$ converges in probability to zero. In contrast, the numerator of $\text{SB}_{\text{BDML}}$ is $\Sigma_{12}$, which has a proper prior that is unaffected by the sample size or the number of control regressors. This is a clear advantage of working with the reduced form: stating a prior in which $\delta$ and $\gamma$ are independent is innocuous because the implied prior for $\delta$ and $\beta = \delta - \alpha \gamma$ are not independent.
Let $\beta^*=\delta^*-\alpha^*\gamma^*$. We have the following result on the posterior mean of $\alpha$ under the na\"ive estimator.
Note that a similar (in)consistency argument can be established for the cases with \(\lambda,p\asymp n\), as (ref) still holds.
\paragraph{Conditional posterior} Note that by Assumption (ref), the true \(\Sigma^*\) is finite and positive definite, so \(\Sigma^{-1}\) is also finite and positive definite in the limit, almost surely in the posterior. Also by Assumption (ref), \(\Sigma_X\) is finite with a positive proportion of positive eigenvalues. If \(X'X\) does not have full rank, its inverse is given by the Moore-Penrose pseudoinverse.
The conditional posterior of \(B\) given \(\Sigma\) is \[ p(\text{vec}(B) \mid X, W, \Sigma) \propto \exp \left(-\frac 1 2 \text{vec}(B - B_n)' V_n^{-1} \text{vec}(B - B_n) \right), \] where
and \(\hat{B} =[\hat\delta,\hat\gamma]=(X'X)^{-1}X'W\) is the OLS estimator of \(B\). And the conditional posterior of \(\Sigma\) given \(B\) is \[ \Sigma \mid X, W, B \sim \text{Inverse-Wishart}\left(\nu_n, \,\Sigma_n\right), \] where
Below, for notation simplicity, let \(B_0=0\).
For the bias-variance trade-off, note that the variance of the posterior mean is not the same as the posterior variance. As we only have posterior variance for the BDML estimator, below we focus on deriving and comparing the variance of the posterior mean.
Under the assumption that $\sigma_{\varepsilon}^{2*}$ and $\sigma_{V}^{2*}$ are known, we can resort to Walker (2024) by verifying his assumptions in our high dimensional regression setup.
Let $P^*$ denote the probability measure based on the true DGP, and $\|\cdot\|_{TV}$ denote the total variation distance.
For a generic matrix $S$, define the Orlicz \(\psi_1\) norm as \[ \|S\|_{\psi_1} = \inf\left\{s>0:\,\mathbb{E}\exp\left(\frac{\|S\|}{s}\right)\le 2\right\}, \] where $\|\cdot\|$ is the spectral norm. Note that $\|S\|_{\psi_1}$ is finite if and only if $S$ is sub-exponential. Then, the matrix Bernstein inequality is as follows.
We now compare our proposed BDML to a number of alternatives from the literature in a simulation experiment adapted from linero_nonparametric_2023. In each replication we generate
We then construct $(D_i, Y_i)$ according to (ref)--(ref) above. This design generates new control regressors $X$ and a new parameter vector $\beta$ in each replication but holds $\gamma$ and $\alpha$ constant across replications. Following the “fixed” design from linero_nonparametric_2023, we vary $\sigma_\varepsilon$ over a grid of values from $1$ to $4$ and set
where $\iota_p$ denotes a $p$-vector of ones.
We consider three groups of estimators. The first group includes two alternative versions of the procedure from (ref), corresponding to alternative prior structures for $(\gamma, \delta)$. Both of our BDML approaches place a flat prior on $\alpha$ and construct a weakly informative prior for $\Sigma$ by placing an LKJ$(4)$ prior on the correlation matrix for $(U,V)$ and independent Cauchy$(0, 2.5)$ priors on both $\sigma_U$ and $\sigma_V$. The first of our BDML approaches, BDML-Basic, places independent $\text{Normal}(0, 5^2)$ on the elements of $\delta$ and $\gamma$. The second, BDML-Hier, allows different standard deviations in the normal shrinkage priors for $\delta$ and $\gamma$. We achieve this via a hierarchical prior that places independent Inverse-Gamma$(2,2)$ hyper-priors on $\sigma_\delta^2$ and $\sigma_\gamma^2$.\footnote{This is equivalent to placing independent Student-t$(4)$ distributions on each of the coefficients $\delta_j$ and $\gamma_j$.} We sample from the posteriors for both BDML approaches via Hamiltonian Monte Carlo (HMC), using the no-U-turn (NUTS) sampler as implemented in the Stan probabilistic programming language.
The second group of estimators we consider contains Bayesian competitors to our BDML approach. This group includes three estimators. The first two, HCPH and Linero, correspond to the proposals from hahn_regularization_2018 and linero_nonparametric_2023 detailed in Algorithms (ref)--(ref) above and rely on a first-stage regression of $D$ on $X$. The third, Na\"ive, estimates a single-equation linear regression model including (ref) only and is a fully-Bayesian variant of the ridge regression approach described in (ref) above. We carry out posterior sampling for these three approaches exactly as in linero_nonparametric_2023. In particular, we rely on the Gibbs sampler implemented in the BLR() function from the R package BLR for Bayesian linear regression. Again following linero_nonparametric_2023 we place a flat prior on $\alpha$ across all three approaches and on $\kappa$ from (ref) and a Jeffreys prior on all error variances. We likewise place default “ridge priors” on the coefficients on $X$ across all regressions in both stages: independent Normal$(0, \sigma^2)$ priors on each coefficient and a Jeffreys prior on $\sigma^2$.\footnote{The resulting prior for the coefficients on $X$ is improper and sharply peaked at zero.}
The third and final group of estimators contains frequentist DML competitors to our BDML approach. Both of the approaches in this group estimate $\alpha$ via (ref) where the “first-step” estimators $\widehat{\gamma}$ and $\widehat{\delta}$ are constructed via ridge regression.\footnote{We use ridge rather than LASSO for these preliminary DML regressions because recent work by kolesar_fragility_2025 and shen_can_2024 suggests that the “double-LASSO” approach belloni_inference_2014 perform poorly due to its fragility in variable selection.} We carry out ridge regression by computing the posterior mean from draws produced using BLR() with the default parameter values, as explained above. We consider two versions of this FDML approach: “full-sample” FDML-Full and “split-sample” FDML-Split. FDML-Full uses the full dataset for both estimation stages whereas FDML-Split randomly splits the sample partitions the sample into two halves, the first of which is used to estimate $(\widehat{\delta}_\text{ML}, \widehat{\gamma}_\text{ML})$ and the second of which is used to residualize $(Y,D)$ and estimate $\alpha$.
For each of the seven estimators of $\alpha$ described above, we compute the root mean-squared error (RMSE), along with the frequentist coverage probability and average width of a corresponding (nominal) 95% confidence interval. For the five Bayesian procedures, interval estimates are constructed as equal-tailed 95% posterior credible intervals. For the two frequentist DML procedures, confidence intervals are constructed based on the usual OLS interval for the regression of the residualized $Y$ on the residualized $D$, following the approach suggested by belloni_inference_2014 and often used in empirical research. Results based on 200 simulation replications from the data-generating process given in (ref)--(ref) appear in (ref) and (ref). Coverage probabilities are very poor for all methods except BDML-Basic, BDML-Hier, and Linero. All three of these approaches produce coverage probabilities close to the nominal 95% level, but BDML-Hier is closest. From panel (ref) of (ref), we also see that BDML-Hier produces the shortest average intervals among the three methods with approximately correct coverage. This holds across all values of $\sigma_\varepsilon$ in our simulation. Turning our attention to RMSE, the three best-performing methods remain BDML-basic, BDML-Hier, and Linero. Once again from panel (ref) of (ref), we see that BDML-Hier performs best. The fact that BDML out-peforms Linero in this setting is particularly encouraging given that this simulation design was taken directly from linero_nonparametric_2023 and favors (ref).
\singlespacing