Bayesian causal inference: a critical review

(a) Basic factorization and versions of causal estimands

Because of the unavoidable missing potential outcomes, causal inference under the potential outcomes framework is inherently a missing data problem [6,15]. The Bayesian paradigm offers a unified framework for statistical inference with missing data and thus for causal inference [27]. Below we review the general structure of Bayesian causal inference that was first outlined in [15].

Four quantities are associated with each unit $i$, $\{Y_i(0),Y_i(1),Z_i,X_i\}$, where $\{Z_i,X_i,Y_i(Z_i)\}$ are observed but $Y_i(1-Z_i)$ is missing. Bayesian inference views all these quantities as random variables and centres around specifying a model for them. Based on the Bayesian model, we can draw inference on causal estimands—functions of the model parameters, covariates and potential outcomes—from the posterior predictive distributions of the parameters and the unobserved potential outcomes. Specifically, we assume the joint distribution of these random variables of all units is governed by a parameter $\theta =({\theta }_X,{\theta }_Z,{\theta }_Y)$, conditional on which the random variables for each unit are i.i.d.. Then we can factorize the joint density $Pr\{Y_i(0),Y_i(1),Z_i,X_i\mid \theta \}$ for each unit $i$ as

Before diving into the technical details, we first clarify the subtle but important difference between the Bayesian estimation of the PATE and SATE estimands. For the PATE, we rewrite the outcome-model-based identification formula in §2 as ${\tau }^{\mathrm{p}}=\int \{{\mu }_1(x;{\theta }_Y)-{\mu }_0(x;{\theta }_Y)\}F(\mathrm{d}x;{\theta }_X)$, which depends only on the unknown parameters ${\theta }_X$ and ${\theta }_Y$. Therefore, Bayesian inference for the PATE requires obtaining posterior distributions of $({\theta }_X,{\theta }_Y)$. By contrast, the SATE ${\tau }^{\mathrm{s}}$ is a function of the potential outcomes ${\{Y_i(0),Y_i(1)\}}_{i=1}^N$, which involves both observed and missing quantities. Bayesian inference for the SATE requires imputing the missing potential outcomes $Y_i^{\text{mis}}$ from their posterior predictive distributions based on the outcome model, and consequently deriving the posterior distribution of ${\tau }^{\mathrm{s}}$.

However, in practice, we rarely model the possibly multi-dimensional covariates $X_i$, and instead condition on the observed values of the covariates. This is equivalent to replacing $F(x;{\theta }_X)$ with ${\hat{\mathbb{F}}}_X$, the empirical distribution of the covariates. Therefore, most Bayesian causal inference (e.g. [9,28]) in fact focuses on the mixed average treatment effect (MATE) [6]

(b) Posterior inference of causal effects

Regardless of the version of the target estimand, the following assumption is commonly adopted.

Assumption 3.2 imposes independent prior distributions for parameters $({\theta }_X,{\theta }_Z,{\theta }_Y)$. It is unique to the Bayesian paradigm of causal inference. It is imposed primarily for modelling and computational convenience and may appear innocuous. However, as elaborated in §4(b), it may lead to unintended and undesirable implications in high-dimensional problems. Under assumptions 2.1 and 3.2, the joint posterior distribution of $\theta =({\theta }_X,{\theta }_Z,{\theta }_Y)$ and the missing potential outcomes is proportional to

By definition, ${\tau }^{\mathrm{p}}=\mathbb{E}\{Y_i(1)\}-\mathbb{E}\{Y_i(0)\}$ does not depend on the association between $Y_i(0)$ and $Y_i(1)$, denoted by the parameter $\rho$. Similarly, $\tau (x)$ does not depend on $\rho$, but ${\tau }^{\mathrm{s}}$ does. So in the inference of ${\tau }^{\mathrm{p}}$ and $\tau (x)$, we usually directly specify the marginal models $Pr\{Y_i(z)\mid X_i;{\theta }_Y\}$ or equivalently $Pr(Y_i\mid Z_i=z,X_i;{\theta }_Y)$ under ignorability [28]. The observed-data likelihood based on (3.4) becomes $\prod _{i:Z_i=1}Pr(Y_i\mid Z_i=1,X_i;{\theta }_Y)\prod _{i:Z_i=0}Pr(Y_i\mid Z_i=0,X_i;{\theta }_Y)$. Imposing a prior for ${\theta }_Y$, we can proceed to infer ${\theta }_Y$ and subsequently ${\tau }^{\mathrm{p}}$, ${\tau }^{\mathrm{m}}$, or $\tau (x)$ using the usual Bayesian inferential procedures.

Bayesian inference of ${\tau }^{\mathrm{s}}$ is more complex, because it depends on both $Y_i(0)$ and $Y_i(1)$ and thus requires posterior sampling of both ${\theta }_Y$ and ${\mathbf{Y}}^{\text{mis}}$. The most common sampling strategy is through data augmentation: iteratively simulate $\theta$ and ${\mathbf{Y}}^{\text{mis}}$ given each other and the observed data, namely from $Pr({\theta }_Y\mid {\mathbf{Y}}^{\text{mis}},{\mathbf{Y}}^{\text{obs}},\mathbf{Z},\mathbf{X})$ and $Pr({\mathbf{Y}}^{\text{mis}}\mid {\mathbf{Y}}^{\text{obs}},\mathbf{Z},\mathbf{X};{\theta }_Y)$. The former, given the observed data and the imputed ${\mathbf{Y}}^{\text{mis}}$, can be straightforwardly obtained by a complete-data analysis based on $Pr\{{\theta }_Y\mid \mathbf{Y}(1),\mathbf{Y}(0),\mathbf{X}\}\propto Pr({\theta }_Y)\prod _{i=1}^{N}Pr\{Y_i(1),Y_i(0)\mid X_i;{\theta }_Y\}$. The latter requires more elaboration. Specifically, we can show that $Pr({\mathbf{Y}}^{\text{mis}}\mid {\mathbf{Y}}^{\text{obs}},\mathbf{Z},\mathbf{X};{\theta }_Y)$ is proportional to $\prod _{i:Z_i=1}Pr\{Y_i(0)\mid Y_i(1),X_i;{\theta }_Y\}\prod _{i:Z_i=0}Pr\{Y_i(1)\mid Y_i(0),X_i;{\theta }_Y\}$. This renders that imputing the ${\mathbf{Y}}^{\text{mis}}$ depends crucially on the joint distribution of $\{Y_i(1),Y_i(0)\}$. Because $Y_i(0)$ and $Y_i(1)$ are never jointly observed, the data provide no information about their association $\rho$. Unless the specific marginal model places constraints on $\rho$, the posterior distribution of $\rho$ would be the same as its prior. Consequently, the posterior distribution of ${\tau }^{\mathrm{s}}$ would be sensitive to the prior of $\rho$.

The above discussion prompts us to clarify the notion of identifiability in Bayesian inference. Under the Frequentist paradigm, a parameter is identifiable if any of its two distinct values give two different distributions of the observed data. Under the Bayesian paradigm, there is no consensus. For example, Lindley [29] argued that all parameters are identifiable in Bayesian analysis because with proper prior distributions, posterior distributions are always proper. In this sense, $\rho$ is identifiable. However, due to the fundamental problem of causal inference, there is no information in the data on $\rho$ and it is reasonable to label it as non-identifiable. This is distinct from the parameters that the data provide direct information on, e.g. those in the marginal distributions of the outcomes in each arm, which are reasonable to label as identifiable. Lindley’s perspective of identifiability blurs such distinction. A more informative perspective is provided by Gustafson [30], who argued that a parameter is weakly or partially identifiable, if a substantial region of its posterior distribution is flat, or its posterior distribution depends crucially on its prior distribution even with large samples, such as $\rho$. Another example of a partially identifiable parameter is $Pr\{Y_i(1)>Y_i(0)\}$, which depends on $\rho$ [6,31]. In this perspective, identifiability in Bayesian inference is no longer all-or-nothing; instead, it is a continuum in between. This issue motivates the strategy of transparent parameterization, where one separates identifiable and non-identifiable parameters, and treats the latter as sensitivity parameters in a sensitivity analysis [32–36]. More discussion will be given in §6(b).

Example 1 revisited We now illustrate the posterior inference of the causal estimands in example 3.1. Here, the parameters $\beta$’s and $\sigma$’s are identifiable, but $\rho$ is not in the Frequentist sense. We fit a Bayesian linear regression model of $Y_i$ on $X_i$ to each observed arm $z$, with independent priors on $({\beta }_1,{\sigma }_1^2)$ and $({\beta }_0,{\sigma }_0^2)$. The observed likelihood factorizes into two parts: the data in treatment group $\{(X_i,Y_i):Z_i=1\}$ and control group $\{(X_i,Y_i):Z_i=0\}$ contribute to the likelihood of $({\beta }_1,{\sigma }_1^2)$ and $({\beta }_0,{\sigma }_0^2)$, respectively. For example, imposing the conventional conjugate normal-inverse ${\chi }^2$ priors, we can draw from the posterior distribution of $\beta$ and $\sigma$, and thus that of the MATE by plugging the posterior draws into the closed-form of ${\tau }^{\mathrm{m}}$ in (3.3). To obtain the PATE, we would have to specify a multivariate model for $Pr(X;\theta )$, and derive the posterior distribution of ${\theta }_X$ and then plug it into the closed form of ${\tau }^{\mathrm{p}}$ in (3.3). This can also be implemented, e.g. via a Bayesian bootstrap step without a model, as described in the next paragraph. To obtain the SATE, we could specify a prior for $\rho$ or fix it to a value. Given $\rho$ and each draw of $({\beta }_1,{\beta }_0,{\sigma }_1^2,{\sigma }_0^2)$, we can impute $Y_i^{\text{mis}}$ as follows: for treated units, $Y_i^{\text{mis}}=Y_i(0)$, and we draw $Y_i(0)$ from $\mathcal{N}({\beta }_0^{\prime }{X}_i+\rho {\sigma }_0/{\sigma }_1\cdot (Y_i-{\beta }_1^{\prime }{X}_i),{\sigma }_0^2(1-{\rho }^2))$; for control units, $Y_i^{\text{mis}}=Y_i(1)$, and we draw $Y_i(1)$ from $\mathcal{N}({\beta }_1^{\prime }{X}_i+\rho {\sigma }_1/{\sigma }_0\cdot (Y_i-{\beta }_0^{\prime }{X}_i),{\sigma }_1^2(1-{\rho }^2)).$ Plugging these posterior predictive draws of $Y_i^{\text{mis}}$ and the observed outcomes into the definition of ${\tau }^{\mathrm{s}}$, we obtain its posterior distribution. We suggest varying the sensitivity parameter $\rho$ from 0 to 1, which corresponds to conditionally independent potential outcomes and perfectly correlated potential outcomes, respectively.

An interesting alternative Bayesian strategy is through the Bayesian bootstrap [37], where the units are re-weighted with weights drawn from a Dirichlet distribution. The Bayesian bootstrap is a general strategy to simulate the posterior distribution of a parameter under a non-parametric model, which can be viewed as the limit of the inference under the Dirichlet Process prior [38]. This renders the Bayesian bootstrap relevant to causal inference in at least two ways. First, one can generate posterior samples from the distribution of $Pr(X_i;{\theta }_X)$ without specifying a parametric model. This is desirable in inferring the population estimands like the PATE and the CATE [39]. However, how to integrate these samples of $X$ into the inference of the target causal estimand is case-dependent and generally adds complexity to the analysis compared to the MATE. Second, the Bayesian bootstrap offers a general recipe for incorporating many standard Frequentist procedures into Bayesian inference. For example, Taddy et al. [40] used it to quantify the uncertainty in linear and tree-based methods for estimating the CATE. Chamberlain & Imbens [41] used it in M-estimation with an application to the setting of instrumental variables (see §7(a)). However, we view the Bayesian bootstrap approach as peripheral to Bayesian causal inference because it does not capitalize on arguably the main strength of Bayesian inference, namely, a unified inferential framework underpinned by the Bayes theorem with versatile choice of priors and outcome models.

(a) Common specification of the outcome model

Section 3 shows that the central component in Bayesian inference of the CATE, PATE, and MATE is to specify the outcome model ${\mu }_z(x)=\mathbb{E}(Y_i\mid X_i=x,Z_i=z;{\theta }_Y)$. We can either model the two treatment groups jointly with a single function ${\mu }_z(x)=\mu (z,x)$ or model each group separately with two functions ${\mu }_1(x)$ and ${\mu }_0(x)$, known as S-learner or T-learner, respectively, in the literature [42]. The most basic outcome model is a linear regression: $\mu (z,x)=x+z+xz$, with Gaussian error terms, where the treatment-covariate interaction term $xz$ captures the treatment effect heterogeneity. This model is equivalent to specifying a linear regression in each group. But the equivalence does not hold for nonlinear models; in fact, S-learners and T-learners with the same type of nonlinear models often lead to markedly different causal estimates.

Linear regressions are easy to implement and interpret, but they are often too restricted. In real-world problems, it is crucial to specify $\mu (z,x)$ flexibly enough to approximate the possibly complex underlying true data generating mechanism. This is particularly desirable as the recent focus in causal inference has been moving toward heterogeneous treatment effects. Outcome modelling is the most natural approach in these studies: one can simply specify an outcome model and derive the CATE as a function of the model parameters. There has been a rapidly increasing adoption of non-parametric and machine learning models for $\mu (z,x)$. One of the most widely used such models is based on regression trees. At a high level, regression trees partition the covariate space into non-overlapping regions and the prediction in each region is based solely on the data that fall in that region. The parameters of a regression tree characterize where to split the covariate space and how to predict the outcomes in a terminal node [43]. An ensemble of regression trees—usually referred to as forests—are often combined to improve the prediction.

Within the Bayesian paradigm, the Bayesian Additive Regression Tree (BART) [44] has become very popular for causal inference. BART places certain priors on the parameters of the regression trees to control the depth of the tree and the degree of shrinkage of the mean parameters in terminal nodes. Hill [9] first advocated using BART to specify the outcome model $\mu (z,x)$ in an S-learner. One can also specify a T-learner with a separate BART model for each treatment group ${\mu }_z(x)$. However, without any additional structure on the marginal models ${\mu }_z(x)$, T-learners often result in large variance of the treatment effects. Hahn et al. [8] proposed the Bayesian Causal Forest method based on an alternative reparametrization, $\mu (z,x)=g_1(x)+g_2(x)z$, where $g_1(x)$ models the distribution of $Y(0)$ and $g_2(x)$ represents the heterogeneous treatment effect, with a separate BART prior for $g_1(x)$ and $g_2(x)$. The BART models have a number of advantages, including fast computation, good performance of default choice of hyperparameters and available software. When a study has adequate covariate overlap, BART has been shown to outperform numerous competing methods, including (the Frequentist) random forests, in many empirical applications, e.g. [45,46]. Other Bayesian non-parametric models, such as Gaussian process [47], Dirichlet process [48–51], have also been considered for causal inference. We refer interested readers to [10] for a more detailed review of these methods.

(b) Challenges in high dimensions

Conducting statistical inference in high dimensions is challenging in general. We differentiate between two high-dimensional settings: (i) an outcome model with an infinite or a large number of parameters, regardless of the number of covariates, such as non-parametric and semiparametric models, and (ii) a large number of covariates. Both settings are increasingly common in causal inference, particularly in studies targeting the CATE. As discussed earlier, outcome modelling is the primary method in these settings, and Bayesian non-parametric models have become a mainstay of the model choice.

A straightforward application of the standard Bayesian non-parametric priors to outcome modelling is sometimes inadequate for causal inference, even with low dimensional covariates. An important consideration is that a desirable prior should accurately reflect uncertainty according to the degree of covariate overlap because intuitively the uncertainty of causal estimates should increase as the degree of overlap decreases. Below, we reproduce a simple example in [52] (first constructed by Surya Tokdar) with a single covariate to illustrate this point.

Example 4.1.

[Choice of priors in estimating the CATE] Consider a study with 250 treated and 250 control units. Each unit has a single covariate $X$ that follows a Gamma distribution with mean 60 and 35 in the control ($Z_i=0$) and treatment ($Z_i=1$) group, respectively, and with s.d. 8 in both groups. To convey the main message, we consider a true outcome model with constant treatment effects: $Y_i(z)=10+5z-0.3X_i+{ϵ}_i$ with ${ϵ}_i\sim \mathcal{N}(0,1)$, where the CATE $\tau (x)=5$ for all $x$. The scatterplots in the upper panel of figure 1 show that covariate overlap is good between the groups in the middle of the range of $X$ (around 40 to 50), but deteriorates towards the tails of $X$.

• Download figure
• Open in new tab
• Download PowerPoint

To estimate the CATE, we fit an outcome model separately in each group: $\mu (z,x)=f_z(x)+{ϵ}_i$ with ${ϵ}_i\sim \mathcal{N}(0,{\sigma }^2)$. We choose three priors for $f_z(x)$: (i) a linear model $f_z(x)={\alpha }_z+{\beta }_{z}x$ with a Gaussian prior for the coefficients; (ii) a BART prior similar to [9]; (iii) a Gaussian Process prior [53] with the covariance function specified by a Gaussian kernel with signal-to-noise ratio parameter $\rho$ and inverse-bandwidth parameter $\lambda$: ${(f_z(x_1),f_z(x_2),\dots ,f_z(x_N))}^{\prime }\sim \mathcal{N}(0,\Sigma )$ where ${\Sigma }_{ij}={\delta }^2{\rho }^2\exp \{-{\lambda }^2||x_i-x_j|{|}^2\})$. Figure 1 shows the posterior means of ${\mu }_z(X)$ and the CATE, with corresponding uncertainty band as a function of $X$. Here, we focus on the uncertainty quantification. In the region of good overlap, all three models lead to similar points and credible interval estimates of the CATE, but a marked difference emerges in the region of poor overlap. The linear model appears overconfident in estimating the CATE. The Gaussian process trades potential bias with wider credible interval as overlap decreases and produces a more adaptive uncertainty quantification. BART produces shorter error bars than the Gaussian Process (but wider than linear models), but the width of the credible interval remains similar regardless of the degree of overlap and thus is over confident in the presence of poor overlap.

Example 4.1 illustrates that, even with low dimensional covariates, standard Bayesian priors can have markedly different operating characteristics when the two groups are poorly overlapped, and not all priors can adaptively capture the uncertainty according to the degree of overlap. A primary reason for BART potentially underestimating uncertainty in poor overlap is its lack of smoothness, in contrast to the Gaussian process. Nonetheless, such a problem can be mitigated by soft decision trees as in [54].

When there are a large number of covariates, the Bayesian paradigm usually achieves regularization through sparsity-inducing priors for the outcome model, such as the spike-and-slab prior [55], the Bayesian LASSO [56], as well as the model averaging techniques [57–59]. The use of these methods in causal inference is surveyed in [10,50]. High-dimensional covariates pose additional complications to Bayesian causal inference. Robins & Ritov [60] pointed out that non-parametric estimates often have slow convergence rates in this regime, which translates into poor finite-sample performance. A central challenge in causal inference is that covariate overlap is rapidly diminishing as the covariates dimension increases, violating the overlap assumption that underpins standard causal analysis [61]. Lack of overlap exacerbates the usual inferential challenges—such as sparsity and slow convergence—in high-dimensional analysis. Even if we assume linear outcome models, we must carefully specify the priors on the regression coefficients. For example, Hahn et al. [7] showed that standard Bayesian regularization on the nuisance parameters may indirectly regularize important causal parameters and thus induce bias, namely the regularization induced confounding. This issue was rigorously investigated in [62]. Specifically, Linero [62] defines the selection bias as ${\delta }_z=\mathbb{E}(Y_i\mid Z_i=z)-\mathbb{E}\{Y_i(z)\}$, and showed that, under the seemingly innocuous prior independence assumption 3.2, many Bayesian regularization priors would a priori induce the selection bias ${\delta }_z$ to sharply concentrate around zero as the number of covariates, $p$, increases, to the extent that no amount of data would overcome such a bias. This implies that assumption 3.2 effectively acts as a strongly informative prior as $p$ increases. Such a phenomenon is referred to as prior dogmaticism and is the Bayesian analogue of the aforementioned problem in Ritov et al. [63]. This line of research highlighted the importance of incorporating the propensity score in Bayesian causal inference [7,62,64,65], which echos the insights from the Frequentist double machine learning method [66,67]. Specifically, the regularized propensity score model or outcome model alone would not be sufficient for valid causal inference, but combining the two would achieve desirable convergence rate and finite sample performance in high-dimensional causal analysis.

(a) Include the propensity score as a covariate in the outcome model

The propensity score was first proposed to be included as the only covariate in a Bayesian outcome model under ignorability, which would reduce the model complexity [68]. However, as later pointed out by Zigler [59]: $Pr\{Y(z)\mid X,e(X)\}=Pr\{Y(z)\mid X\}\ne Pr\{Y(z)\mid e(X)\}$. So using the propensity score as the single covariate in the outcome model would not lead to the target outcome distribution $Pr\{Y(z)\mid X\}$, but using it as an additional covariate, i.e. specifying a model $\mu (z,x)=\mu (z,x,e(x))$, would. This specification is effectively conducting an outcome regression at each value of the propensity score, and thus can be viewed as a smoothed version of combining propensity score stratification and outcome modelling. In a sense, this specification provides a Bayesian analogue of the double robustness [52,69]. On the one hand, when the outcome model is correctly specified, $\mu (Z_i,X_i,e(X_i))$ reduces to $\mu (Z_i,X_i)$ because $e(X_i)$ is a function of $X_i$ and thus is redundant regardless of its specification. On the other hand, when the outcome model is misspecified but the propensity model is correctly specified, the results are robust to the outcome model specification because the treatment and control groups are approximately balanced in covariates within each stratum of the propensity score. Various reparametrizations have been proposed. One example is to specify $\mu (z,x,e(x))=g_1(x,e(x))+g_2(z,x)$, with $g_1(\cdot )$ being a non-parametric model and $g_2(\cdot )$ being a parametric model. Little [70] adopted a penalized spline model of $e(x)$ for $g_1(\cdot )$. In the aforementioned Bayesian Causal Forest, Hahn et al. [8] imposed a separate BART model for $g_1(\cdot )$ and $g_2(\cdot )$, and demonstrated that adding the propensity score as an additional predictor in $g_1$ significantly improves the empirical estimation of the CATE.

This strategy is usually implemented in two stages: first estimate the propensity score as $\hat{e}(X)$ and then plug it into the Bayesian outcome model $\mu (Z,X,\hat{e}(X))$. Such a two-stage procedure is not dogmatically Bayesian, which generically refers to the procedure of specifying a model with parameters and prior distributions of these parameters and then use the Bayes theorem to obtain the posterior distributions of the parameters. A direct consequence is that this procedure may not properly propagate the uncertainty of estimating the propensity score in the outcome model [69]. A dogmatic Bayesian approach would jointly model $e(X;{\theta }_Z)$ and $\mu (Z,X,e(X);{\theta }_Y)$ and draw posterior inference of ${\theta }_Z$ and ${\theta }_Y$ simultaneously [71]. However, when the outcome model is misspecified, the joint-modelling approach would introduce a feedback problem, that is, the fit of the outcome model would inform the estimation of the propensity scores. This violates the unconfoundedness assumption, distorts the balancing property of the propensity score, and consequently leads to biased estimate of causal effects. A suggested remedy is to first fit a Bayesian model for $e$ and then plug the posterior predictive draws of $e$ into the outcome [11]. Such a two-stage procedure is still not dogmatically Bayesian, but provides more robust posterior inference to model misspecification empirically.

However, adding the propensity score into the outcome model is controversial conceptually, because the outcome model reflects the nature of the generating process of the potential outcomes, which arguably should not depend on how the treatment is assigned [72].

(b) Dependent priors

The Bayesian causal inference outlined in §3 rests on the prior independence assumption 3.2, without which the propensity score model cannot be ignored from the likelihood. But this assumption is not always plausible in real applications. Various priors that do not rely on this assumption have been proposed [47,63–65]. Below we show two simple examples.

The first example is due to [58] and is designed for simultaneous variable selection for the propensity score and outcome models. Specifically, assume a logistic propensity score model $\text{logit}\{Pr(Z_i=1\mid X_i)\}={\alpha }^{\prime }{X}_i$ and a linear outcome model $Y_i\mid Z_i,X_i\sim \mathcal{N}(\tau Z_i+{\beta }^{\prime }{X}_i,{\sigma }^2)$. Assume each of the $j$th components of the coefficients, ${\alpha }_j$, follow the spike-and-slab prior [73]: ${\alpha }_j\mid {\gamma }_j^{\alpha }\sim (1-{\gamma }_j^{\alpha })I_0+{\gamma }_j^{\alpha }\mathcal{N}(0,{\sigma }_{\alpha }^2)$, where ${\gamma }_j^{\alpha }$ is a latent indicator of whether $X_j$ is included in the model and $I_0$ denotes the point mass at 0. A similar spike-and-slab prior is assumed for the coefficients of the outcome model with a latent inclusion indicator ${\gamma }_j^{\beta }$: ${\beta }_j\mid {\gamma }_j^{\beta }\sim (1-{\gamma }_j^{\beta })I_0+{\gamma }_j^{\beta }\mathcal{N}(0,{\sigma }_{\beta }^2)$. Then assume the probability of the events $\{{\gamma }_j^{\alpha }=0\}$ and $\{{\gamma }_j^{\beta }=0\}$ are dependent a priori: $Pr({\gamma }_j^{\beta }=1\mid {\gamma }_j^{\alpha }=1)/Pr({\gamma }_j^{\beta }=0\mid {\gamma }_j^{\alpha }=1)=\omega$, where $\omega \in [1,\mathrm{\infty })$ is a dependence hyperparameter that controls the prior odds of including $X_j$ into the outcome model when it is included in the propensity score model. Larger $\omega$ implies stronger prior dependence between the variable selection in the two models.

The second example is due to [74]. Assume $Y_i(1)\mid X_i\sim \mathcal{N}({\mu }_1,{\sigma }_1^{2}e(X_i))$ and $Y_i(0)\mid X_i\sim \mathcal{N}({\mu }_0,{\sigma }_0^2(1-e(X_i)))$, with flat priors on ${\mu }_1$ and ${\mu }_0$. If the propensity scores are known, the posterior mean of the PATE equals the Hajék version of the IPW estimator:

Carefully designed dependent priors often achieve desirable finite sample results and are more reasonable in real world studies. However, specification of such priors is case-dependent, and there is no general solution.

(c) Posterior predictive inference

A general, albeit not dogmatically Bayesian, strategy is to specify both a propensity score model $e(X_i;{\theta }_Z)$ and an outcome model $\{{\mu }_1(X_i;{\theta }_Y),{\mu }_0(X_i;{\theta }_Y)\}$, and obtain posterior draws of $e(X_i;{\theta }_Z)$ and $\{{\mu }_1(X_i;{\theta }_Y),{\mu }_0(X_i;{\theta }_Y)\}$ from their respective posterior predictive distributions, and then plug these posterior draws into the doubly-robust estimator ${\hat{\tau }}^{\text{dr}}$ [66,75]. A variance estimator of the resulting estimator ${\hat{\tau }}^{\text{dr}}$ is given in [66]. In the same vein, Ding & Guo [76] incorporated the propensity score in Bayesian posterior predictive $p$-value. For the model with the Fisher’s sharp null hypothesis of no causal effect for any units whatsoever (i.e. $Y_i(1)=Y_i(0)$ for all $i$), the procedure in [76] is equivalent to the Fisher randomization test averaged over the posterior predictive distribution of the propensity score. Simulations in [76] show the advantages of the Bayesian $p$-value compared to the Frequentist analogue. This perspective offers a straightforward and flexible strategy to integrate Bayesian modelling and common Frequentist procedures for causal inference and enables proper uncertainty quantification.

Besides the above three strategies, another general approach is through the aforementioned Bayesian bootstrap, which can be used to simulate the posterior distribution of any parameter that can be formulated as M-estimation or estimating equation [41,77]. As special cases, because the IPW estimator ${\hat{\tau }}^{\text{ipw}}$ and the doubly robust estimator ${\hat{\tau }}^{\text{dr}}$—both involving the propensity scores—are both solutions to estimating equations, they can be naturally combined with the Bayesian bootstrap to devise a Bayesian version. However, such an approach may be guilty of ‘Bayesian for the sake of being Bayesian’, and their methodological and practical value compared to competing methods is unclear.

(a) Parametrization involving distributions with unmeasured confounders

The first parametrization used for sensitivity analysis is motivated by the intuition that a hidden confounder may completely explain away the association between the treatment and the outcome even after adjusting for observed covariates. In a historic debate, Fisher [78] hypothesized that the strong association between cigarette smoking and lung cancer might be due to a hidden genetic factor as their ‘common cause’ or confounder. Cornfield et al. [79] derived an inequality showing that to explain away the observed association, the association between the unmeasured confounder and cigarette smoking must be larger than or equal to the association between cigarette smoking and lung cancer. Their work helped to initiate the field of sensitivity analysis.

Let $U$ denote an unmeasured confounder and assume that unconfoundedness holds conditional on $(X,U)$: $Z\perp \perp \{Y(0),Y(1)\}\mid X,U$. The joint distribution of all variables factorizes into

As an extension of [79], Ding & VanderWeele [84] treated the treatment-confounder ($Z$ and $U$) and outcome-confounder ($Y$ and $U$) associations as two sensitivity parameters, and derived analytical thresholds for them in order to explain away the observed treatment-outcome ($Z$ and $Y$) association. Based on that theory, VanderWeele & Ding [85] further simplified by assuming the two associations to be the same and called the resulting threshold the E-value, as a measure of robustness of the causal conclusions with respect to unmeasured confounding. The E-value framework is model-free because it avoids modelling assumptions with $U$; it also avoids repeating the analysis over a range of sensitivity parameters as in the competing methods and thus is simple to calculate.

(b) Parametrization involving distributions of potential outcomes

The second parametrization is motivated by an alternative mathematical expression of the unconfoundedness assumption: $Pr\{Y(z)\mid Z=1,X\}=Pr\{Y(z)\mid Z=0,X\}$ for $z=0,1$, representing the fact that the units in the two randomized arms are comparable in terms of potential outcomes. This class of sensitivity analysis is based on sensitivity parameters that directly represent the difference between the distributions $Pr\{Y(z)\mid Z=1,X\}$ and $Pr\{Y(z)\mid Z=0,X\}$ instead of modelling the difference with an unobserved $U$. This is implemented in the context of time-varying treatments (see §7(b)) and Frequentist semiparametric estimation [86]. Franks et al. [34] pointed out the importance of distinguishing between model fit and sensitivity to unconfoundedness: the former involves identifiable parameters (e.g. the parameters in the model of the marginal distributions of the outcome $Pr\{Y_i(z)\mid Z_i=z,X\}$) whereas the latter involves unidentifiable parameters (e.g. the association between $Y_i(1)$ and $Y_i(0)$). The merit of this parametrization is apparent in this perspective because it separates identifiable and unidentifiable parameters. Franks et al. [34] proceeded under the Bayesian paradigm and used a copula—parameters of which are the sensitivity parameters—to connect the two identifiable marginal distributions.

A related branch of sensitivity analysis is Rosenbaum’s bounds [87]. His original formulation takes the association between $Z$ and the potential outcomes conditional on observed $X$, denoted by $\Gamma$, as the sole sensitivity parameter for quantifying unmeasured confounding. He has also made connections to the parametrization in §6(a) [88]. Starting with a matched sample to mimic a randomized matched-pairs experiment, one can then repeat the Fisher randomization test on the sharp null hypothesis of no treatment effect given a range of $\Gamma$ values, and find the threshold of $\Gamma$ at which the $p$-value of the test changes from significant to insignificant. This approach was later generalized to derive the bounds of a given estimator under different $\Gamma$ values. Grounded in Fisherian randomization inference, Rosenbaum’s framework does not have a natural Bayesian analogue.

Besides the above two classes, there are numerous other approaches to sensitivity analysis based on alternative parameterization of unmeasured confounding. However, a common criticism of various approaches to sensitivity analysis is that, in order to assess the consequence of the untestable unconfoundedness assumption, one has to make even more untestable assumptions, e.g. specifying models involving $U$. Moreover, sensitivity analysis, after all, is a secondary analysis in causal studies, and thus simple implementation and intuitive interpretation is much desired. The considerations underpin the dominance of the E-value method over other methods in practice, particularly in medicine and public health.

(a) From instrumental variable to principal stratification

Instrumental variable (IV) is one of the most important techniques for causal inference in economics and social sciences. IVs are used in settings where dependence of the assignment on the potential outcomes cannot plausibly be ruled out, even conditional on observed covariates. An IV is a variable that provides a source of exogenous (or unconfounded) variation that helps identify causal effects. IV methods are based on a set of assumptions alternative to ignorability. Specifically, an IV satisfies three conditions: (i) it occurs before a treatment; (ii) it is independent of the treatment-outcome confounding; and (iii) it affects the outcome only through its (non-zero) effects on the treatment assignment. Finding a valid IV is challenging in observational studies and many clever natural experiments have been identified [89]. Given a valid IV, one can extract the causal effects of the treatment on an outcome by a two-stage least-squares (2SLS) estimator: first, fit a linear regression of the treatment on the IV; second, fit a linear regression of the outcome on the fitted value of the treatment from the first stage, the coefficient of which is taken as the causal effect of the treatment on the outcome. Covariates can be added in both stages.

The IV method has been developed within the structural equation model framework (see [89] for a review), and the 2SLS IV estimator may not correspond to a causal effect within the potential outcomes framework except for a few special cases. In a landmark paper, Angrist et al. [90] connects IV to the potential outcomes framework in the setting of randomized experiments with binary treatment and all-or-nothing compliance, with the initial random assignment playing the role of an IV. But many questions remain on the connection between the IV method and the potential outcomes framework in more general settings. Below we will describe the special setting of Angrist et al. [90].

We introduce some new notation. For unit $i$, let $Z_i$ be the randomly assigned treatment (1 for the treatment and 0 for the control), and $W_i$ be the actual treatment status (1 for the treatment and 0 for the control). When $Z_i\ne W_i$, non-compliance arises. Because $W_i$ occurs post-assignment, it has two potential values, $W_i(0)$ and $W_i(1)$, with $W_i=W_i(Z_i)$. As before, the outcome $Y_i$ has two potential outcomes, $Y_i(0)$ and $Y_i(1)$. Based on their joint potential status of the actual treatment $U_i=(W_i(1),W_i(0))$, the units fall into four compliance types: compliers $U_i=(1,0)=\text{co}$, never-takers $U_i=(0,0)=\text{nt}$, always-takers $U_i=(1,1)=\text{at}$ and defiers $U_i=(0,1)=\text{df}$ [90]. A key property of $U_i$ is that it is not affected by the treatment assignment, and thus can be regarded as a pre-treatment characteristic. Therefore, comparisons of $Y_i(1)$ and $Y_i(0)$ within the stratum of $U_i$ have standard subgroup causal interpretations: ${\tau }_u=\mathbb{E}\{Y_i(1)-Y_i(0)\mid U_i=u\}$, for $u=\text{nt},\text{co},\text{at},\text{df}$; ${\tau }_u$ are later called principal causal effects. The conventional causal estimand in clinical trials is the intention-to-treat effect that ignores the compliance information, which is the weighted average of the four stratum-specific effects: $\mathbb{E}\{Y_i(1)-Y_i(0)\}=\sum _u{\pi }_u{\tau }_u$, where ${\pi }_u=Pr(U_i=u)$ is the proportion of the stratum $u$. The intention-to-treat effect measures the effect of the assignment instead of the actual treatment.

Due to the fundamental problem of causal inference, individual compliance stratum $U_i$ is not observed. So the principal causal effects are non-identifiable without additional assumptions. Besides randomization of $Z_i$, Angrist et al. [90] make two additional assumptions: (i) monotonicity: $W_i(1)\ge W_i(0)$, and (ii) exclusion restriction: $Y_i(1)=Y_i(0)$ whenever $W_i(1)=W_i(0)$. Monotonicity rules out defiers, and exclusion restriction imposes that the assignment has zero effects among never-takers and always-takers. Then the compiler average causal effect is identified by

We now describe the Bayesian inference of the IV set-up, first outlined in [91]. Without additional assumptions, the observed cells of $Z$ and $W$ consist of a mixture of units from more than one stratum. For example, the units who are assigned to the treatment arm and took the treatment ($Z=1,W=1$) can be either always-takers or compliers. One must disentangle the causal effects for different compliance types from observed data. Therefore, model-based inference here resembles that of a mixture model. In Bayesian analysis, it is natural to impute the missing label $U_i$ under some model assumptions. Specifically, six quantities are now associated with each unit, $\{Y_i(1)$, $Y_i(0)$, $W_i(1)$, $W_i(0)$, $X_i$, $Z_i\}$, four of which, $\{Y_i^{\text{obs}}=Y_i=Y_i(Z_i),W_i^{\text{obs}}=W_i=W_i(Z_i),Z_i,X_i\}$, are observed and the remaining two, $\{Y_i^{\text{mis}}=Y_i(1-Z_i),W_i^{\text{mis}}=W_i(1-Z_i)\}$, are unobserved. Assume the joint distribution of these random variables of all units is governed by a parameter $\theta$, conditional on which the random variables for each unit are iid. We assume unconfoundedness $Pr\{Z_i=1\mid X_i,W_i(1),W_i(0),Y_i(1),Y_i(0)\}=Pr(Z_i=1\mid X_i)$, and impose a prior distribution $Pr(\theta )$. Then the joint posterior distribution of $\theta$ and the missing potential outcomes are proportional to the complete-data likelihood as follows:

Using the same arguments as in §3, to infer population and mixed estimands, we only need to specify two marginal outcome models for $Y_i(1)$ and $Y_i(0)$ instead of a joint model, and do not need to impute the missing potential outcomes $Y^{\text{mis}}$. But we do need to impute the latent $U_i$, or, equivalently, the missing intermediate variable $W^{\text{mis}}$. We can simulate the joint posterior distribution $Pr(\theta ,{\mathbf{W}}^{\text{mis}}\mid {\mathbf{Y}}^{\text{obs}},{\mathbf{W}}^{\text{obs}},\mathbf{Z},\mathbf{X})$ by iteratively imputing the missing $\mathbf{U}$ from $Pr({\mathbf{W}}^{\text{mis}}\mid {\mathbf{Y}}^{\text{obs}},{\mathbf{W}}^{\text{obs}},\mathbf{Z},\mathbf{X},\theta )$ and updating the posterior distribution of $\theta$ from $Pr(\theta \mid {\mathbf{Y}}^{\text{obs}},{\mathbf{W}}^{\text{obs}},{\mathbf{W}}^{\text{mis}},\mathbf{Z},\mathbf{X})$.

Below we illustrate the Bayesian procedure via a simple example of the IV approach.

Frangakis & Rubin [93] generalized the IV approach to principal stratification, a unified framework for causal inference with post-treatment confounding. In the simplest scenario, a post-treatment confounded variable lies in the causal pathway between the treatment and the outcome; it cannot be adjusted in the same fashion as a pre-treatment covariate in causal inference. A principal stratification with respect to a post-treatment variable is the classification of units based on the joint potential values of the post-treatment variable, and the stratum-specific effects are called principal causal effects, of which ${\tau }_{\text{co}}$ is a special case. The post-treatment variable setting includes a wide range of examples. For instance, in the non-compliance setting, the ‘treatment’ is the randomized treatment assignment, the ‘post-treatment’ variable is the actual treatment received, and the compliance types are the principal strata [91,92,94,95]. Zeng et al. [96] connects principal stratification to the local IV method with a continuous IV and binary treatment [97]. Other examples include censoring due to death [98], surrogate endpoints [99,100], regression discontinuity designs [101], time-varying treatments [102], and many more. The choices of target strata and thus estimands, interpretations, and identifying assumptions depend on specific applications, details of which are omitted here.

(b) Time-varying treatment and confounding

In real-world situations, subjects often receive treatments sequentially at multiple time points, and the treatment assignment at each time is affected by both baseline and time-varying confounders as well as previous treatments [103,104]. Such settings are referred to as time-varying, or sequential, or longitudinal treatments.

Consider a study where treatments are assigned at $T$ time points. Let $Z_{it}$ denote the treatment at time $t$ for unit $i$ $(i=1,\dots ,N;t=1,\dots ,T)$. At baseline ($t=0$), each unit $i$ has time-invariant covariates $L_{i0}$ measured; then after the treatment assignment at time $t-1$ and prior to the assignment at time $t$, a set of time-varying confounders $L_{i,t-1}$ are measured, which include the intermediate measurements of the final outcome and the covariates that are affected by the previous treatments. For example, in a cancer study, baseline covariates can include sex, age, race and time-varying confounders can include intermediate cancer progression and other clinical traits such as blood pressure measured prior to the next treatment. Denote the observed and hypothetical treatment sequence of length $t$ by ${\overline{Z}}_{it}=(Z_{i1},\dots ,Z_{it})$ and ${\overline{z}}_t=(z_1,\dots ,z_t)$, respectively, and the sequence of time-varying confounders by ${\overline{L}}_{it}=(L_{i0},L_{i1},\dots ,L_{it})$. For each ${\overline{z}}_T$, there is a potential outcome $Y({\overline{z}}_T)$. The final observed outcome $Y_i=Y_i({\overline{Z}}_{iT})$, corresponding to the entire observed treatment sequence, is measured after treatment assigned at $T$. A common causal estimand is the marginal effect comparing two pre-specified treatment sequences, $\overline{z},{\overline{z}}^{\prime }\in {\{0,1\}}^T$: ${\tau }_{{\overline{z}}_T,{\overline{z}}_T^{\prime }}=\mathbb{E}\{Y_i({\overline{z}}_T)-Y_i({\overline{z}}_T^{\prime })\}$. For simplicity, below we drop the subscript $i$.

The central question to causal inference with sequential treatments is the role of the time-varying confounders $L_t$ in the assignment mechanism. These variables are affected by the previous treatments and also affect the future treatment assignment and outcome. Much of the literature assumes a sequentially ignorable assignment mechanism [103], that is, the treatment at each time is unconfounded conditional on the observed history, which consists of past treatments ${\overline{Z}}_{t-1}$ and time-varying confounders ${\overline{L}}_{t-1}$, as stated below.

A full Bayesian approach to time-varying treatments [105] would specify a joint model for treatment assignment $Z_t$ and time-varying confounders $L_t$ at all time points as well as all the potential outcomes $Y({\overline{z}}_T)$, and then derive the posterior predictive distributions of the missing potential outcomes and thus of the estimands. This procedure is a straightforward extension from the structure introduced in §3. However, the joint modelling approach is rarely used because it quickly becomes intractable as the time $T$ and the number of time-varying confounders increases.

Instead, most of the Bayesian methods are grounded in the g-computation. Under sequential ignorability, the average potential outcome $\mathbb{E}\{Y_i({\overline{z}}_T)\}$ is identified from the observed data via the $g$-formula [103]: