A Bayesian view of doubly robust causal inference

TL;DR

This paper provides a Bayesian perspective on doubly robust causal inference, highlighting limitations of likelihood-based methods and proposing a posterior predictive approach that decouples outcome and treatment models for improved robustness.

Contribution

It introduces a Bayesian posterior predictive method that incorporates inverse treatment probabilities as importance sampling weights, offering a new way to achieve doubly robust causal inference.

Findings

Simulations demonstrate the proposed method's effectiveness.

Likelihood-based Bayesian inferences cannot fully explain doubly robust properties.

Decoupling models enhances robustness in causal effect estimation.

Abstract

In causal inference confounding may be controlled either through regression adjustment in an outcome model, or through propensity score adjustment or inverse probability of treatment weighting, or both. The latter approaches, which are based on modelling of the treatment assignment mechanism and their doubly robust extensions have been difficult to motivate using formal Bayesian arguments, in principle, for likelihood-based inferences, the treatment assignment model can play no part in inferences concerning the expected outcomes if the models are assumed to be correctly specified. On the other hand, forcing dependency between the outcome and treatment assignment models by allowing the former to be misspecified results in loss of the balancing property of the propensity scores and the loss of any double robustness. In this paper, we explain in the framework of misspecified models why doubly robust inferences cannot arise from purely likelihood-based arguments, and demonstrate this through simulations. As an alternative to Bayesian propensity score analysis, we propose a Bayesian posterior predictive approach for constructing doubly robust estimation procedures. Our approach appropriately decouples the outcome and treatment assignment models by incorporating the inverse treatment assignment probabilities in Bayesian causal inferences as importance sampling weights in Monte Carlo integration.

Tables

Table 1. Table 1: Estimates for the marginal causal contrast (true value = 1 absent 1 =1 ) over 1000 simulation rounds

Scenario	Estimator	Point	Relative	SD	SE	MC	Coverage
		estimate	bias (%)			error	(%)
(I)	Naive	0.347	$-$65.3	0.128	0.128	0.004	0.3
	Adjusted	0.667	$-$33.3	0.091	0.092	0.003	3.6
	IPTW	1.001	0.1	0.135	0.138	0.005	94.4
	OR/PS (obs. information)	0.997	$-$0.3	0.071	0.095	0.002	99.3
	OR/PS (adj. sandwich)	0.997	$-$0.3	0.071	0.073	0.002	95.3
	DR	0.998	$-$0.2	0.087	0.088	0.003	94.2
	Clever covariate	1.026	2.6	0.110	0.110	0.003	93.0
	OR/IPTW	0.990	$-$1.0	0.082	0.083	0.003	93.9
	Two-step (forward sampling)	0.999	$-$0.1	0.071	0.113	0.002	99.8
	Two-step (variance decomposition)	0.995	$-$0.5	0.071	0.112	0.002	99.8
	Joint estimation	1.046	4.6	0.071	0.071	0.002	89.5
	Importance sampling	0.991	$-$0.9	0.083	0.080	0.003	93.6
	Importance sampling/DR	0.997	$-$0.3	0.087	0.086	0.003	93.9
(II)	Naive	0.347	$-$65.3	0.128	0.128	0.004	0.3
	Adjusted	0.997	$-$0.3	0.065	0.067	0.002	95.3
	IPTW	0.629	$-$37.1	0.133	0.131	0.005	20.3
	OR/PS (obs. information)	0.997	$-$0.3	0.071	0.071	0.002	95.4
	OR/PS (adj. sandwich)	0.997	$-$0.3	0.071	0.071	0.002	95.4
	DR	0.999	$-$0.1	0.074	0.074	0.003	95.9
	Clever covariate	0.999	$-$0.1	0.075	0.075	0.003	95.9
	OR/IPTW	0.999	$-$0.1	0.074	0.073	0.003	95.7
	Two-step (forward sampling)	1.000	0.0	0.071	0.072	0.002	96.5
	Two-step (variance decomposition)	0.997	$-$0.3	0.071	0.071	0.002	95.5
	Joint estimation	0.997	$-$0.3	0.071	0.071	0.002	95.4
	Importance sampling	0.999	$-$0.1	0.074	0.072	0.003	95.2
	Importance sampling/DR	0.999	$-$0.1	0.074	0.073	0.003	95.5

Equations

$$\mathrm{E}({\mathbf{Y}}_{1i})-\mathrm{E}({\mathbf{Y}}_{0i})={\mathrm{E}}_{X_i}\{\mathrm{E}({\mathbf{Y}}_{1i}|X_i)\}-{\mathrm{E}}_{X_i}\{\mathrm{E}({\mathbf{Y}}_{0i}|X_i)\}.$$

$$\mathrm{E}({\mathbf{Y}}_{1i})-\mathrm{E}({\mathbf{Y}}_{0i})={\int }_{x_i}\left\{\mathrm{E}(Y_i|Z_i=1,x_i)-\mathrm{E}(Y_i|Z_i=0,x_i)\right\}p(x_i)\text{d}{x}_i.$$

$$p(v)$$

$${\int }_{\varphi }\frac{1}{n}\sum _{i=1}^n\left\{m(1,x_i;\varphi )-m(0,x_i;\varphi )\right\}{\pi }_n(\varphi )\text{d}\varphi ,$$

$$\frac{1}{n}\sum _{i=1}^n\left[m\left\{1,e(b_i;\hat{\gamma }),s_i;\hat{\varphi }\right\}-m\left\{0,e(b_i;\hat{\gamma }),s_i;\hat{\varphi }\right\}\right],$$

$$\frac{1}{n}\sum _{i=1}^n\frac{y_i{z}_i-\{z_i-e(b_i;\hat{\gamma })\}m(1,s_i;\hat{\varphi })}{e(b_i;\hat{\gamma })}$$

$$-\frac{1}{n}\sum _{i=1}^n\frac{y_i(1-z_i)-[(1-z_i)-\{1-e(b_i;\hat{\gamma })\}]m(0,s_i;\hat{\varphi })}{1-e(b_i;\hat{\gamma })},$$

$$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\bigg{\{}\frac{z_{i}}{e(b_{i};\widehat{\gamma})}-\frac{1-z_{i}}{1-e(b_{i};\widehat{\gamma})}\bigg{\}}\left\{y_{i}-m(z_{i},s_{i};\widehat{\phi})\right\}+\frac{1}{n}\sum_{i=1}^{n}\bigg{\{}m(1,s_{i};\widehat{\phi})-m(0,s_{i};\widehat{\phi})\bigg{\}}.$$

$$m(z,s;\varphi )={\varphi }_0+{\varphi }_{1}z+{\varphi }_2^{\top }s+{\varphi }_{3}c(z,b).$$

$$\frac{1}{n}\sum _{i=1}^n\left\{m(1,x_i;\hat{\varphi })-m(0,x_i;\hat{\varphi })\right\}={\hat{\varphi }}_1+{\hat{\varphi }}_3\frac{1}{n}\sum _{i=1}^n\left\{\frac{1}{e(b_i;\hat{\gamma })}-\frac{1}{1-e(b_i;\hat{\gamma })}\right\}.$$

$$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\bigg{\{}$$

$$\frac{1}{n}\sum _{i=1}^n\left\{\mathrm{E}({\mathbf{Y}}_{1i}|s_i;\hat{\varphi })-\mathrm{E}({\mathbf{Y}}_{0i}|s_i;\hat{\varphi })\right\},$$

$$l_i(\varphi )=\sum _{a=0}^1{1}_{\{z_i=a\}}\frac{\log p({\mathbf{y}}_{ai}|s_i;\varphi )}{\mathrm{pr}(Z_i=a|b_i)}.$$

$$u_i(\varphi )$$

$${\int }_{\psi ,\gamma ,\varphi }{\int }_{x_i}m(a,e(b_i;\gamma ),s_i;\varphi )p(x_i;\psi ){\pi }_n(\varphi |\gamma ){\pi }_n(\gamma ){\pi }_n(\psi )\text{d}{x}_i\text{d}\varphi \text{d}\gamma \text{d}\psi ,$$

$${\pi }_n(\varphi |\gamma )\propto \prod _{i=1}^{n}p\{y_i|z_i,e(b_i;\gamma ),s_i;\varphi \}{\pi }_0(\varphi )\text{ and }{\pi }_n(\gamma )\propto \prod _{i=1}^{n}p(z_i|b_i;\gamma ){\pi }_0(\gamma ).$$

$$\mathrm{var}({\Phi }_1|v)$$

$$p(v)={\int }_{{\varphi }^{\ast },{\gamma }^{\ast },\psi }\left\{\prod _{i=1}^{n}p(y_i,z_i|x_i;{\varphi }^{\ast },{\gamma }^{\ast })p(x_i;\psi )\right\}{\pi }_0({\varphi }^{\ast }){\pi }_0({\gamma }^{\ast }){\pi }_0(\psi )\text{d}{\varphi }^{\ast }\text{d}{\gamma }^{\ast }\text{d}\psi ,$$

$${\pi }_n({\varphi }^{\ast },{\gamma }^{\ast })\propto \prod _{i=1}^n\left[p\{y_i|z_i,e(b_i;{\gamma }^{\ast }),s_i;{\varphi }^{\ast }\}p(z_i|b_i;{\gamma }^{\ast })\right]{\pi }_0({\varphi }^{\ast }){\pi }_0({\gamma }^{\ast }).$$

$${\int }_{{\varphi }^{\ast },{\gamma }^{\ast },\psi }{\int }_{x_i}m\{a,e(b_i;{\gamma }^{\ast }),s_i;{\varphi }^{\ast }\}p(x_i;\psi ){\pi }_n({\varphi }^{\ast },{\gamma }^{\ast }){\pi }_n(\psi )\text{d}{x}_i\text{d}{\varphi }^{\ast }\text{d}{\gamma }^{\ast }\text{d}\psi .$$

$${\mathrm{E}}_{\mathcal{E}}\{l(\varphi ;V_i)|v\}={\mathrm{E}}_{\Xi |V}\left[{\mathrm{E}}_{\mathcal{E}}\{l(\varphi ;V_i)|v;\xi \}\right],$$

$${\mathrm{E}}_{\mathcal{E}}[l(\varphi ;V_i)|v;\xi ]$$

$$={\int }_{v_i}l(\varphi ;v_i)\frac{p_{\mathcal{E}}(v_i|v;\xi )}{p_{\mathcal{O}}(v_i|v;\xi )}{p}_{\mathcal{O}}(v_i|v;\xi )\text{d}{v}_i$$

$$={\int }_{v_i}l(\varphi ;v_i)\frac{p_{\mathcal{E}}(y_i|z_i,x_i,v;\xi )p_{\mathcal{E}}(z_i)p_{\mathcal{E}}(x_i|v;\xi )}{p_{\mathcal{O}}(y_i|z_i,x_i,v;\xi )p_{\mathcal{O}}(z_i|x_i,v;\xi )p_{\mathcal{O}}(x_i|v;\xi )}{p}_{\mathcal{O}}(v_i|v;\xi )\text{d}{v}_i$$

$$={\int }_{v_i}l(\varphi ;v_i)\frac{p_{\mathcal{E}}(z_i)}{p_{\mathcal{O}}(z_i|x_i,v;\xi )}{p}_{\mathcal{O}}(v_i|v;\xi )\text{d}{v}_i,$$

$${\mathrm{E}}_{\mathcal{E}}[l(\varphi ;V_i)|v;\xi ]={\int }_{v_i}l(\varphi ;v_i)w_i(\xi )\sum _{k=1}^n{\xi }_k{\delta }_{v_k}(v_i)\text{d}{v}_i=\sum _{k=1}^n{w}_k(\xi ){\xi }_{k}l(\varphi ;v_k),$$

$$\hat{\varphi }(\xi )\equiv \arg \underset{\varphi }{\max }\sum _{k=1}^n{w}_k(\xi ){\xi }_{k}l(\varphi ;v_k).$$

$${\mathrm{E}}_{\Xi |V}\{{\mathrm{E}}_{\mathcal{E}}(Y_i|Z_i=a,v;\xi )\}$$

$$={\int }_{\xi }\sum _{k=1}^n{\xi }_{k}m\{a,s_k;\hat{\varphi }(\xi )\}{p}_{\mathcal{O}}(\xi |v)\text{d}\xi ,$$

$${\mathrm{E}}_{\mathcal{E}}(Y_i|Z_i=1,v;\xi )-{\mathrm{E}}_{\mathcal{E}}(Y_i|Z_i=0,v;\xi )$$

$$=\sum _{k=1}^n{\xi }_k\{y_i-m(z_k,x_k;\xi )\}\left\{\frac{z_k}{{\mathrm{pr}}_{\mathcal{O}}(Z_k=1|x_k,v;\xi )}-\frac{1-z_k}{{\mathrm{pr}}_{\mathcal{O}}(Z_k=0|x_k,v;\xi )}\right\}$$

$$+\sum _{k=1}^n{\xi }_k\left\{m(1,x_k;\xi )-m(0,x_k;\xi )\right\},$$

$${\mathrm{E}}_{\mathcal{E}}(Y_i|Z_i=a,v)$$

$$=\frac{{\int }_{\varphi ,\psi }{\int }_{x_i}{\int }_{y_i}{y}_{i}p(y_i|Z_i=a,x_i;\varphi ){\mathrm{pr}}_{\mathcal{E}}(Z_i=a)p(x_i;\psi ){\pi }_n(\varphi ){\pi }_n(\psi )\text{d}{y}_i\text{d}{x}_i\text{d}\varphi \text{d}\psi }{{\int }_{\varphi ,\psi }{\int }_{x_i}{\int }_{y_i}p(y_i|Z_i=a,x_i;\varphi ){\mathrm{pr}}_{\mathcal{E}}(Z_i=a)p(x_i;\psi ){\pi }_n(\varphi ){\pi }_n(\psi )\text{d}{y}_i\text{d}{x}_i\text{d}\varphi \text{d}\psi }$$

Full text

A revised version of this article has been accepted for publication in Biometrika, published by Oxford University Press.

Saarela, O., Belzile, L. R. and D. A. Stephens. A Bayesian view of doubly robust causal inference, Biometrika (2016), 103 (3): 667-681. doi:10.1093/biomet/asw025.

1. Introduction

In causal inference contexts, confounding is most often controlled by one of two approaches: first, by conditioning on the confounders in a regression model for the outcome in an outcome regression model; secondly, by modelling the treatment assignment mechanism to obtain so-called propensity score values, and then using these scores to construct strata within the observed sample, or a pseudo-sample from a hypothetical population, within which treatment assignment is not confounded. This pseudo-sample can be obtained via inverse probability of treatment weighting of the original sample, analogously to survey sampling procedures. The outcome regression adjustment method requires correct specification of the regression function in order to obtain consistent inference; this may be achieved in practice using flexible regression strategies and complex functions of typically a large number of covariates. The propensity score adjustment methods focus principally on the specification of the treatment assignment model, which may be similarly flexible or complex. Under either approach, sufficient control of confounding is therefore dependent on possibly unverifiable modelling assumptions. This has motivated the development of doubly robust methods in which both model components are specified, but only one of them needs to be correctly specified to sufficiently control for confounding.

Adjustments that depend on the propensity score using regression or reweighting are not easy to interpret from the Bayesian perspective, since Bayesian inferences are naturally based on modelling of the outcome, with modelling of the treatment assignment playing no role in inference relating to the outcome/treatment relationship (Robins and Ritov, , 1997). Gustafson, (2012) attempted a Bayesian interpretation as a compromise between a saturated outcome model and a parametric one; however, the treatment assignment model did not feature in this interpretation. Scharfstein et al., (1999) and Bang and Robins, (2005) pointed out a connection between a doubly robust estimator and a model-based estimator; the two are equivalent if the outcome model in the latter features the so called clever covariate, a particular function of the propensity score.

Separately from the above developments, there has been a body of work studying Bayesian versions of propensity score adjustment to control for confounding (Hoshino, , 2008; McCandless et al., 2009a, ; McCandless et al., 2009b, ; McCandless et al., , 2010, 2012; An, , 2010; Kaplan and Chen, , 2012; Zigler et al., , 2013; Chen and Kaplan, , 2015; Graham et al., , 2015). These approaches require one of two kinds of compromises; they either force a parametrization which makes the outcome and treatment assignment models dependent, thus losing the balancing property of the propensity score, or cut such a feedback in which case the inference procedures are no longer Bayesian. Because of these difficulties, some authors (e.g. Achy-Brou et al., , 2010, p. 828) have been content to fix the propensity scores to their best estimates in model-based inferences, without attempting to jointly estimate the two model components. In an alternative approach, Wang et al., (2012) and Zigler and Dominici, (2014) have suggested connecting the outcome and treatment assignment models through the prior distribution in order to incorporate the uncertainty in confounder selection. Herein we do not consider model uncertainty, but rather, concentrate on inferences with a priori specified outcome and treatment assignment models.

The purpose of this paper is to clarify the theoretical and practical motivations for Bayesian propensity score adjustment, and the relationships between the different methods proposed for this, which have not been fully explored previously. We address these in the context of double adjustment for both the potential confounders and the propensity score, and argue that the problem cannot be properly understood without considering it in the framework of misspecified models. To provide an alternative to Bayesian propensity score adjustment, we propose deriving Bayesian versions of various inverse probability of treatment weighted estimators, including inverse probability of treatment weighted outcome regression and the semi-parametric double robust estimator, through posterior predictive expectations, with the weights introduced as importance sampling weights in Monte Carlo integration.

2. Preliminaries

2. Notation and assumptions

For simplicity, we consider the case of a single binary treatment, and defer discussion of the longitudinal, multiple exposure and continuous cases to the Discussion. Let the random vectors $X_{i}$ represent a set of pre-treatment covariate measurements, $Z_{i}$ a binary treatment allocation indicator, and $Y_{i}$ an outcome for individual $i$, measured after sufficient time has passed since administering the treatment. We adopt, for convenience, the standard potential outcome (or counterfactual) framework: for individual $i$, the observed outcome is related to the two possible potential outcomes $(\mathbf{Y}_{0i},\mathbf{Y}_{1i})$ by $Y_{i}=(1-Z_{i})\mathbf{Y}_{0i}+Z_{i}\mathbf{Y}_{1i}$. We assume that $X_{i}$ includes a sufficient set of covariates to control for confounding in the sense that Z_{i}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}(\mathbf{Y}_{0i},\mathbf{Y}_{1i})\,|\,X_{i} (cf. ignorable treatment assignment, Rosenbaum and Rubin, , 1983, p. 43). The propensity score $e(X_{i})\equiv\mathrm{pr}(Z_{i}=1\,|\,X_{i})$ has the balancing property Z_{i}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X_{i}\,|\,e(X_{i}), which also implies that (\mathbf{Y}_{0i},\mathbf{Y}_{1i})\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z_{i}\,|\,e(X_{i}) – this scalar score is therefore useful in controlling for confounding.

If the covariate space $X_{i}$ is high-dimensional, in practice the task of controlling for confounding often involves some covariate selection; here one can either select for the features that predict the outcome, or the treatment assignment. To represent this, let $S_{i}$ and $B_{i}$ represent some a priori selected subsets of the all observed features $X_{i}$, so that the latter can be partitioned as $X_{i}=(S_{i},R_{i})$ or $X_{i}=(B_{i},C_{i})$, where possibly $R_{i}=\emptyset$ and $C_{i}=\emptyset$. If the selected set of features $S_{i}$ captures all relevant prognostic information, then \mathbf{Y}_{0i}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X_{i}\,|\,S_{i} (Hansen, , 2008). For the remainder of the paper, we consider the stronger condition (i) (\mathbf{Y}_{0i},\mathbf{Y}_{1i})\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X_{i}\,|\,S_{i}, which requires that $S_{i}$ also captures all relevant information about possible effect modification.

In this case $S_{i}$ is sufficient to control for confounding, since from the properties of conditional independence (Dawid, , 1979) it follows that (\mathbf{Y}_{0i},\mathbf{Y}_{1i})\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z_{i}\,|\,S_{i}. If, on the other hand, the selected set of features $B_{i}$ has the balancing property (ii) Z_{i}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X_{i}\,|\,B_{i}, it is sufficient to control for confounding. We are interested in estimation procedures which are valid when either (\mathbf{Y}_{0i},\mathbf{Y}_{1i})\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X_{i}\,|\,S_{i} (but possibly $Z_{i}{\not\perp\!\!\!\!\perp}X_{i}\,|\,B_{i}$) or Z_{i}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X_{i}\,|\,B_{i} (but possibly $(\mathbf{Y}_{0i},\mathbf{Y}_{1i}){\not\perp\!\!\!\!\perp}X_{i}\,|\,S_{i}$).

Our parameter of interest is an average causal contrast such as

[TABLE]

Interest then lies in the identifiability of the average potential outcomes based on the observed data. When the ‘no unmeasured confounder’ assumption holds, it follows that (e.g. Hernán and Robins, , 2006, p. 43)

[TABLE]

2. Bayesian model formulation for outcome regression

Any Bayesian model specification is constructed via de Finetti’s representation for exchangeable random sequences $v_{i}\equiv(x_{i},y_{i},z_{i})$ (see for example Bernardo and Smith, , 1994, Chap. 4). The (subjective) joint distribution for observations $v\equiv(v_{1},\ldots,v_{n})$ may then be represented by

[TABLE]

implying the existence of the parametric models and the prior density $\pi_{0}(\phi,\gamma,\psi)$. Since the representation theorem is not constructive, that is, does not specify the models implicit in (1), in practice inferences about a given finite-dimensional parametrization involves the often implicit assumption that $p(y_{i}\,|\,z_{i},x_{i};\phi_{0})=f(y_{i}\,|\,z_{i},x_{i})$, $p(z_{i}\,|\,x_{i};\gamma_{0})=f(z_{i}\,|\,x_{i})$ and $p(x_{i};\psi_{0})=f(x_{i})$, where $(\phi_{0},\gamma_{0},\psi_{0})$ is the limiting value of the posterior $\pi_{n}(\phi,\gamma,\psi)\equiv p(\phi,\gamma,\psi\,|\,v)$ in the sense of e.g. van der Vaart, (1998, p. 139), and where the $f$s represent the ‘true’ limiting (sampling) distributions. We might further assume that the parameters are a priori independent, so that $\pi_{0}(\phi,\gamma,\psi)=\pi_{0}(\phi)\pi_{0}(\gamma)\pi_{0}(\psi)$, in which case it it also follows also that the posterior factorizes as $\pi_{n}(\phi,\gamma,\psi)=\pi_{n}(\phi)\pi_{n}(\gamma)\pi_{n}(\psi)$ (e.g. Gelman et al., , 2004, p. 354–355).

The marginal distribution $p(x_{i};\psi)$ can in practice be specified nonparametrically and estimated using the empirical covariate distribution, leading to a Bayesian estimator of the average causal contrast,

[TABLE]

where $\pi_{n}(\phi)\propto\textstyle\prod_{i=1}^{n}p(y_{i}\,|\,z_{i},x_{i};\phi)\pi_{0}(\phi){\,\textrm{d}}\phi$, and $m(z,x;\phi)\equiv\int yp(y\,|\,z,x;\phi){\,\textrm{d}}y$. In Supplementary Appendix 1, we show that (2) can be motivated without the use of potential outcomes notation through posterior predictive expectations for a new observation under a hypothetical completely randomized setting.

The estimator (2) is the Bayesian version of the well-known direct standardization or g-formula. We return to it in Section 6, but note here that it does not feature the treatment assignment model; rather, the posterior predictive approach for estimating the marginal causal contrast depends entirely on correct specification of the distribution $p(y_{i}\,|\,z_{i},x_{i};\phi)$, or in a moment-based representation, $m(z,x;\phi)$. Before discussing Bayesian alternatives to (2), we briefly review some commonly used frequentist approaches for combining outcome regression and propensity score adjustment.

3. Frequentist approaches for combining outcome regression and propensity score

adjustment

3. Including propensity scores into outcome regression

Because of the balancing property of the propensity score, it is tempting to specify a propensity score $e(b_{i};\gamma)\equiv\mathrm{pr}(Z_{i}=1\,|\,b_{i};\gamma)$ and use a statistical model such as $p\{y_{i}\,|\,z_{i},e(b_{i}),s_{i};\phi\}$, in the hope that, if the prognostic model is is misspecified, adjusting for the propensity score would still sufficiently control for any residual confounding. For simplicity, we take the parameters $\phi$ to specify also the functional dependence between the propensity score and the outcome; to model this dependency, it is advisable to use flexible formulations such as splines (e.g. Zhang and Little, , 2009). Using such an outcome model, the marginal causal contrast would then be estimated by

[TABLE]

where $m\{z,e(b;\gamma),s;\phi\}\equiv\int yp\{y\,|\,z,e(b;\gamma),s;\phi\}{\,\textrm{d}}y$, and where $\widehat{\phi}$ and $\widehat{\gamma}$ are maximum likelihood estimators for the parameters in the outcome regression and treatment assignment model, respectively. The motivation for such a double adjustment is that it is sufficient to control for confounding if either condition (i) or (ii) of Section 2.1 applies. We summarize this property in the following theorem (a proof in Supplementary Appendix 2; in the notation all convergencies are in probability unless otherwise stated).

Theorem 1.

*Estimator (3) is consistent if the outcome model is correctly specified in the sense that $m\{z,e(b;\gamma_{0}),s;\phi_{0}\}=\int yf(y\,|\,z,e(b),s){\,\textrm{d}}y$, the parameters in this can be consistently estimated so that $\widehat{\phi}\rightarrow\phi_{0}$, the treatment assignment model is correctly specified in the sense that $p(z_{i}\,|\,b_{i};\gamma_{0})=f(z_{i}\,|\,b_{i})$ and $\widehat{\gamma}\rightarrow\gamma_{0}$, and either (i) holds true, or (ii) holds true. *

The estimator (3) may be considered ‘doubly robust’ in terms of the covariate selection in the sense that only one of the sets $S_{i}$ and $B_{i}$ needs to be correctly specified, although it still always relies on a correct parametric specification of the model for the expected outcome, conditional on $\{Z_{i},e(B_{i}),S_{i}\}$.

3. The clever covariate and augmented outcome regression

The estimator discussed in the previous section did not specify which function of the propensity score $e(B_{i})$ should be added to the regression model. Scharfstein et al., (1999, p. 1141–1142) and Bang and Robins, (2005, p. 964–965) drew a connection between propensity score regression adjustment and doubly robust estimators of the form

[TABLE]

which can be equivalently represented as

[TABLE]

On considering the score equation derived from a regression of $Y_{i}$ on $Z_{i}$ and $S_{i}$ with mean function $m(z,s;\phi)$, this form suggests incorporating the derived covariate $c(z_{i},b_{i})=z_{i}/e(b_{i})-(1-z_{i})/\{1-e(b_{i})\}$ (termed the ‘clever covariate’ by Rose and van der Laan, , 2008, p. 8) additively into the outcome regression, that is, for example

[TABLE]

The first term in (4) is then zero through the maximum likelihood score equation, leaving only the last term which is the model based estimator of the marginal treatment effect. Thus with the clever covariate in the outcome model, the doubly robust estimator is equivalent to the model-based estimator. In the special case of model (5), this becomes

[TABLE]

A potential drawback of using this covariate is that it may lead to extreme variability for the resulting mean difference estimator, even compared to inverse probability of treatment weighted estimators of the form

[TABLE]

To see why this is the case, the distribution of $c(z_{i},b_{i})$ in itself can be very skewed to the right, but this becomes even more pronounced in the model-based estimator where the clever covariate has to evaluated at both $c(1,b_{i})$ and $c(0,b_{i})$ for each $i=1,\ldots,n$. In contrast, the inverse probability of treatment weighted estimator (6) involves only the probabilities of treatments that were actually assigned.

Eq. (4) is doubly robust in a stronger, semi-parametric, sense than (3); it does not require correct specification of the outcome model, if the treatment assignment model is correctly specified. The approximate Bayesian double robust approach proposed by Graham et al., (2015) involved replacing $m(z,x;\phi)$ in (2) with a linear predictor augmented with the clever covariate. We take this to be a special case of the two-step Bayesian methods to be discussed in Section 4, and thus do not separately consider it in the present paper. However, in Section 6.2 we will show how the form (4) may be derived through posterior predictive expectations and importance sampling.

3. Inverse probability of treatment weighted outcome regression

Yet another estimator for the marginal causal contrast is

[TABLE]

where the parameters $\phi$ in the model for the potential outcomes $(\mathbf{Y}_{1i},\mathbf{Y}_{0i})$ are estimated using an inverse probability of treatment weighted estimating function $l(\phi)=\sum_{i=1}^{n}l_{i}(\phi)$, where

[TABLE]

The corresponding estimating equation is $u(\phi)=\sum_{i=1}^{n}u_{i}(\phi)=0$, where the pseudo-score function is

[TABLE]

Here the treatment assignment probabilities $\mathrm{pr}(Z_{i}=a\,|\,b_{i})$ would in practice be replaced with estimates $\mathrm{pr}(Z_{i}=a\,|\,b_{i};\widehat{\gamma})$. To motivate the use of (8) for the parameter estimation, we present the following theorem (a proof given in Supplementary Appendix 2).

Theorem 2.

*The estimating equation $u(\phi)=0$ is unbiased if the model for the potential outcomes is correctly specified in the sense that $p(\mathbf{y}_{ai}\,|\,s_{i};\phi_{0})=f(\mathbf{y}_{ai}\,|\,s_{i})$, and either (i) or (ii) holds true. *

If we further assume the consistency of the estimator for $\phi$, as well as consistency of $\widehat{\gamma}$ when the weights are correctly specified, it also follows that (7) consistently estimates the marginal causal contrast. In Section 6.1 we demonstrate that an estimator of the form (7) can also be motivated from Bayesian arguments.

4. Two-step estimation when the propensity score is unknown

In observational settings the function $e(B_{i})$ is unknown and has to be estimated. When using an estimator of the form (3), a central question from a Bayesian perspective then is how the uncertainty in the estimation of the parameters $\gamma$ is incorporated in the inference of the marginal causal contrast. A ‘Bayesian’ approach could be motivated by writing the posterior predictive expectation as

[TABLE]

where

[TABLE]

The integrals of the form (9) could be evaluated by Monte Carlo integration, by forward sampling first from $\pi_{n}(\gamma)$ and given the current value $\gamma$, from the conditional posterior $\pi_{n}(\phi\,|\,\gamma)$. However, a concern related to such an approach is that the product of the posterior distributions $\pi_{n}(\phi\,|\,\gamma)$ and $\pi_{n}(\gamma)$ in (9) does not necessarily correspond to any well defined joint posterior, except in the case when the outcome is correctly specified in the sense (i), in which case it does not depend on the propensity score. In contrast, for the correctly specified models in (1) we indeed have that $p(\phi\,|\,v)p(\phi\,|\,x,z)=p(\phi,\gamma\,|\,v)$. As a result, the above outlined two-step approach is not proper Bayesian, and would have to be evaluated on its frequency-based properties.

To give an example of a situation where the two-step Bayesian approach does not result in correct frequency-based inferences, we can consider the special case of the outcome model $m\{z,e(b;\gamma),s;\phi\}=\phi_{0}+\phi_{1}z+\phi_{2}^{\top}s+\phi_{3}^{\top}g\{e(b;\gamma)\}$, where $g$ is, for example, some appropriate spline basis transformation of the propensity score. With this model the estimator based on (9) for the average causal contrast $\mathrm{E}(\mathbf{Y}_{1i})-\mathrm{E}(\mathbf{Y}_{0i})$ reduces to $\int_{\phi,\gamma}\phi_{1}\pi_{n}(\phi\,|\,\gamma)\pi_{n}(\gamma){\,\textrm{d}}\phi{\,\textrm{d}}\gamma$, that is, to an estimator of the posterior mean $\mathrm{E}_{\Gamma\,|\,X,Z}\left\{\mathrm{E}(\Phi_{1}\,|\,v;\gamma)\right\}$. This estimator in turn can be approximated with $\sum_{j=1}^{m}\widehat{\phi}_{1}(\gamma^{(j)})/m$, where $(\gamma^{(1)},\ldots,\gamma^{(m)})$ is a Monte Carlo sample from $\pi_{n}(\gamma)$ (cf. Kaplan and Chen, , 2012, p. 589). We note first that this estimator has the same asymptotic distribution as $\widehat{\phi}_{1}(\widehat{\gamma})$, where the treatment assignment model parameters have been fixed to their maximum likelihood estimates (see Supplementary Appendix 3). In Supplementary Appendix 3 we further show that $\mathrm{avar}\{\widehat{\phi}_{1}(\widehat{\gamma})\}\leq\mathrm{avar}\{\widehat{\phi}_{1}(\gamma_{0})\}$, where $\widehat{\phi}_{1}(\gamma_{0})$ is the estimator given the ‘true’ propensity scores. Thus, with the propensity score adjusted outcome model specification, a variance adjustment due to estimating the propensity scores should reduce the asymptotic variance of the resulting treatment effect estimator compared to a hypothetical situation where the true propensity scores are known (Henmi and Eguchi, , 2004, cf.). In contrast, Kaplan and Chen, (2012, p. 592) and Graham et al., (2015, p. 11) propose variance estimators based on the variance decomposition formula

[TABLE]

which appears to add a further variance component. An explanation for the discrepancy is that with the correctly specified models in the representation (1), we have that $p(\phi_{1}\,|\,v;\gamma)=p(\phi_{1}\,|\,v)$, and the the second variance component becomes zero. On the other hand, if the models are misspecified, the product form likelihood in the representation (1) does not apply in the first place. This illustrates the difficulty in applying Bayesian procedures in the context of incompatible models. Even though this is routinely done in the context of multiple imputation (e.g. Rubin, , 1996), and often produces reasonable results, in the present context there is little motivation to use an approach which introduces an additional variance component to the posterior variance, given that estimation of the propensity scores should reduce the variance of the treatment effect estimator. We further discuss this discrepancy in the following section.

5. Joint estimation of outcome and treatment assignment models

Even though the outcome $Y_{i}$ can obviously be predictive of the individual treatment assignment $Z_{i}$, the outcomes are not informative of the treatment assignment mechanism (a proof in Supplementary Appendix 2):

Theorem 3.

*If the outcome model is correctly specified, the outcomes are non-informative of the parameters characterizing the treatment assignment process. *

In such a case the treatment assignment model plays no part in the inferences, since the corresponding posterior predictive estimator is (2). However, the Bayesian propensity score approach proposed by McCandless et al., 2009a specifies a parametrization making the outcome and treatment assignment models dependent and estimates the parameters jointly. More recently, Zigler et al., (2013) suggested that a similar approach could be used to obtain a Bayesian analogue to doubly robust inferences. Such an approach can be understood by assuming that there exists a de Finetti parametrization $(\phi^{*},\gamma^{*})$ for which

[TABLE]

where $p(y_{i},z_{i}\,|\,x_{i};\phi^{*},\gamma^{*})=p\{y_{i}\,|\,z_{i},e(b_{i};\gamma^{*}),s_{i};\phi^{*}\}p(z_{i}\,|\,b_{i};\gamma^{*})$. Compared to $(\phi,\gamma)$ in (1), neither $\phi^{*}$ or $\gamma^{*}$ retains the original interpretation, but now there is a well defined joint posterior distribution

[TABLE]

Inferences could now be based on the posterior predictive expectations

[TABLE]

At first sight, (11) would seem more natural than (9), since the specification (11) does not make use of incompatible models. However, now the quantities $e(b_{i};\gamma^{*})$ do not possess the balancing properties of propensity scores, and thus it would be difficult to show whether (11) would have the ‘double robustness’ property of the estimator (3).

To address the lack of balance, McCandless et al., (2010) suggested a Gibbs sampler type approach similar to that of Lunn et al., (2009) to cut the feedback from the outcome model, by successively drawing from the conditional posteriors $\pi_{n}(\gamma)$ and $\pi_{n}(\phi\,|\,\gamma)$ to approximate the joint posterior of $(\phi^{*},\gamma^{*})$. However, as discussed in the previous Section, these posteriors are incompatible and generally such a sampling procedure is not guaranteed to converge to any well defined joint distribution. In fact, if the conditional posteriors can be sampled directly, or if the second sampling step is allowed to converge to the corresponding conditional distribution, the inferences based on the formulations (9) and (11) will be equivalent.

To sum up, trying to recover fully probabilistic inferences through sampling from a joint posterior of the outcome and treatment assignment model parameters loses the balancing property of the propensity scores, and consequently, the properties of the resulting estimator would be difficult to establish. On the other hand, cutting the feedback in an attempt to recover the balancing property would mean that the inferences are no longer based on well-defined posterior distributions. Thus, in the following section, following the approach outlined in Saarela et al., 2015b , we formulate alternative Bayesian estimators that are not based on Bayesian propensity score adjustment.

6. Posterior predictive inferences with importance sampling

6. Inverse probability of treatment weighted outcome regression

It has been recognized by various authors (e.g. Røysland, , 2011; Chakraborty and Moodie, , 2013) that inverse probability of treatment weighting can be motivated through a change of probability measures, or equivalently, importance sampling. However, as far as we know, before Saarela et al., 2015b this approach has not been used to formulate Bayesian causal inferences. Here we argue that this approach can be used to resolve the paradoxes discussed in Sections (4) and (5). We follow the posterior predictive reasoning of Supplementary Appendix 1, but rather than trying to directly predict a new observation under the experimental setting $\mathcal{E}$, we consider first the task of parameter estimation under this regime. For this purpose, a Bayes estimator for the outcome model parameters can be constructed by maximizing the expected utility $\mathrm{E}_{\mathcal{E}}\{l(\phi;V_{i})\,|\,v\}$ with respect to $\phi$, where the log-likelihood $l(\phi;v_{i})\equiv\log p(y_{i}\,|\,z_{i},s_{i};\phi)$ takes the role of a parametric utility function, and the expectation is over a predicted new observation $v_{i}=(y_{i},z_{i},x_{i})$, $i\notin\{1,\ldots,n\}$.

Let further $\xi$ be a set of parameters characterizing the entire data-generating mechanism under the observational regime $\mathcal{O}$. We can further write the expectation as

[TABLE]

where, following Walker, (2010, p. 26–27), we can consider the lower dimensional decision of maximizing the expected utility $\mathrm{E}_{\mathcal{E}}\{l(\phi;V_{i})\,|\,v;\xi\}$ with respect to $\phi$ conditional on $\xi$. With a known regime $\mathcal{E}$ and the stability assumption discussed in Supplementary Appendix 1, $\arg\max_{\phi}\mathrm{E}_{\mathcal{E}}\{l(\phi;V_{i})\,|\,v;\xi\}$ is a deterministic function of $\xi$. Thus, the uncertainty represented by the posterior distribution $p_{\mathcal{O}}(\xi\,|\,v)$ then also reflects the uncertainty on $\phi$, providing means to construct a posterior distribution for $\phi$. This proceeds as follows; the inner expectation can be written as

[TABLE]

where in the last equality we made use of the stability assumption of Supplementary Appendix 1. In the last form we can replace the predictive distribution under $\mathcal{O}$ with the Bayesian bootstrap specification $p_{\mathcal{O}}(v_{i}\,|\,v;\xi)=\textstyle\sum_{k=1}^{n}\xi_{k}\delta_{v_{k}}(v_{i})$, where $\xi\equiv(\xi_{1},\ldots,\xi_{n})$ and $\Xi\,|\,v\sim\textrm{Dirichlet}(1,\ldots,1)$, as in the weighted likelihood bootstrap of Newton and Raftery, , 1994. Denoting $w_{i}(\xi)\equiv p_{\mathcal{E}}(z_{i})/p_{\mathcal{O}}(z_{i}\,|\,x_{i},v;\xi)$, the expected utility now becomes

[TABLE]

that is, a weighted log-likelihood, motivating the estimator

[TABLE]

An approximate posterior distribution for $\phi$ under $\mathcal{E}$ can now be constructed by repeatedly sampling the weight vectors from $\Xi\,|\,v\sim\textrm{Dirichlet}(1,\ldots,1)$, and recalculating $\widehat{\phi}(\xi)$ at each realization. This approach of creating a mapping between the non-parametric specification and a parametrization relevant to inferences is analogous to Chamberlain and Imbens, (2003), but adds the importance sampling weights to the Dirichlet weights in order to make inferences across the observational and experimental regimes.

The weights $w_{i}(\xi)$ function as importance sampling weights in Monte Carlo integration; they add variability to the estimation, but in the present context provide some protection towards misspecification of the outcome model, in the sense of Section 3.3. The above did not yet address how to calculate these weights; in principle these are fully determined by the current realization of $\xi$ under the non-parametric specification, but in practice parametric model specifications are needed for smoothing purposes, and we need a way to link the $\xi$ and the treatment assignment model parameters $\gamma$. For this purpose $\gamma$ itself can be estimated through the weighted likelihood bootstrap since this readily gives the deterministic function linking the two parametrizations; thus in (12) we choose $w_{i}(\xi)=p_{\mathcal{E}}(z_{i})/p_{\mathcal{O}}\{z_{i}\,|\,b_{i};\widehat{\gamma}(\xi)\}$, where $\widehat{\gamma}(\xi)\equiv\arg\max_{\gamma}\sum_{k=1}^{n}\xi_{k}\log p(z_{k}\,|\,b_{k};\gamma)$. The probabilities $p_{\mathcal{E}}(z_{i})$ are given by the chosen regime $\mathcal{E}$ that is the object of inference; in practice the estimation is most efficient when we choose the target regime to be as close as possible to the observed regime $\mathcal{O}$; this can be achieved by fixing $p_{\mathcal{E}}(z_{i})$ to the marginal treatment assignment probabilities under $\mathcal{O}$, which would result in the usual kind of stabilized inverse probability of treatment weights used in marginal structural modelling (Robins et al., , 2000; Hernán et al., , 2001; Cole and Hernán, , 2008).

The marginal causal contrast may now be estimated through the expectations

[TABLE]

where we used the non-parametric specification $p(s_{i}\,|\,v;\xi)=\sum_{k=1}^{n}\xi_{k}\delta_{s_{k}}(s_{i})$, and where again $p_{\mathcal{O}}(\xi\,|\,v)$ is replaced with the uniform Dirichlet distribution. Since all uncertainty is contained in the parameter vector $\xi$, a posterior distribution for the predictive mean or mean difference can be constructed as before through resampling. The point estimator given by (13) is the direct Bayesian analogue of (7), where the outcome model was estimated using inverse probability of treatment weighted regression. In fact, if we fix $\xi_{k}=1/n$, $k=1,\ldots,n$, instead of considering these as unknown parameters, the two estimators are equivalent. Thus, we conjecture that the estimator given by (13) has a similar ‘double robustness’ property as (7). We demonstrate this through simulations in Section 7, but before that, we show how the semi-parametric doubly robust estimator (4) can be motivated as a posterior predictive expectation.

6. Doubly robust estimation

In Supplementary Appendix 4 we show that under the non-parametric specification in terms of $\xi$, the posterior predictive causal contrast may be written as

[TABLE]

which corresponds to formulation (4). Since the non-parametric specification places no restrictions on the conditional distributions, in practice, to obtain an estimator, the non-parametrically specified quantities $m(z_{k},x_{k};\xi)$ and $\mathrm{pr}_{\mathcal{O}}(Z_{k}=a\,|\,x_{k},v;\xi)$ would have to be replaced with the parametric versions $m\{z_{k},s_{k};\widehat{\phi}(\xi)\}$ and $\mathrm{pr}_{\mathcal{O}}\{Z_{k}=a\,|\,b_{k};\widehat{\gamma}(\xi)\}$, estimated using the weighted likelihood bootstrap, as in the previous section. It is straightforward to see that if the outcome model is correctly specified, the expression (6.2) reduces to the model-based estimator (the second additive term). This again reflects the fact that if we believe in our outcome model, the treatment assignment model does not play a part in the inferences. However, a Bayesian might want to use an estimator of the form (6.2) if being restricted by two parametric constraints, in terms of $\phi$ and $\gamma$, but not knowing which one of these is correct. If either one of the parametric specifications is correct, (6.2) will still reduce to the posterior predictive mean difference, the natural Bayesian estimator. We summarize this property in the following theorem (a proof in Supplementary Appendix 2).

Theorem 4.

*The estimator obtained by substituting the parametric specifications $m\{z_{k},s_{k};\widehat{\phi}(\xi)\}$ and $\mathrm{pr}_{\mathcal{O}}\{Z_{k}=a\,|\,b_{k};\widehat{\gamma}(\xi)\}$ into expression (6.2) is equivalent to the posterior predictive mean difference if either one of these models is correctly specified. *

A posterior distribution for the mean difference can be generated as in the previous section through resampling of the parameter vectors $\xi$ and recalculating (6.2) at each realization; we will illustrate this in the following section.

7. Simulation study

Above we have made a distinction between model misspecification due to omission of relevant covariates, and misspecification of the parametric functional relationship between the outcome and the covariates, and noted that all the estimators discussed in Section 3 should be ‘doubly robust’ against the former type of misspecification. However, in practice the consequences of these two types of misspecification will often be similar; they result in residual confounding. Therefore, herein we investigate how the different estimators perform when the covariate sets $S_{i}$ and $B_{i}$ are not only created by a partitioning, but also a transformation of the $x_{i}$’s. For this purpose, we simulated $X_{ij}\sim\mathrm{N}(0,1)$, $j=1,\ldots,4$, and transformed these as $c_{ij}=|x_{ij}|/({1-2/\pi})^{1/2}$. The true treatment assignment and outcome mechanisms were specified as $Z_{i}\,|\,x_{i}\sim\mathrm{Bernoulli}(\textrm{expit}\{0.4c_{i1}+0.4x_{i2}+0.8x_{i3}\})$ and $Y_{i}\,|\,z_{i},x_{i}\sim\mathrm{N}(z_{i}-c_{i1}-x_{i2}-x_{i4},1)$, respectively. For estimation, we considered two scenarios: (I) $s_{i}\equiv(x_{i1},x_{i2},x_{i3})$ and $b_{i}\equiv(c_{i1},x_{i2},x_{i4})$ (misspecified outcome model and correctly specified treatment assignment model), and (II) $s_{i}\equiv(c_{i1},x_{i2},x_{i3})$ and $b_{i}\equiv(x_{i1},x_{i2},x_{i4})$ (correctly specified outcome model and misspecified treatment assignment model).

We are interested in the marginal causal contrast $\mathrm{E}(\mathbf{Y}_{i1})-\mathrm{E}(\mathbf{Y}_{i0})=1$. To estimate this, we applied the Bayesian estimators discussed in Sections 4, 5, and 6. In the two-step estimation we attempted both forward sampling from the posterior distributions, and the variance decomposition formula (10). In the former, instead of Markov chain Monte Carlo, we applied normal approximations for the posterior distributions, of the form $\Phi\,|\,v;\gamma\sim\mathrm{N}\{\widehat{\phi}(\gamma),S\}$, where $\widehat{\phi}(\gamma)$ is the maximum likelihood estimate and $S$ its estimated variance-covariance matrix. The posterior distribution $\Gamma\,|\,x,z$ was approximated using the weighted likelihood bootstrap. In the joint estimation, we again used a normal approximation, centered at the joint maximum likelihood estimates $(\widehat{\phi},\widehat{\gamma})$, and the variance-covariance matrix given by the inverse of the observed information at the maximum likelihood point. In both two-step and joint estimation, the fitted models were specified as $m\{z_{i},e(b_{i};\gamma),s_{i};\phi\}=\phi_{0}+\phi_{1}z_{i}+\phi_{2}^{\top}s_{i}+\phi_{3}^{\top}g\{e(b_{i};\gamma)\}$, where $g$ is a cubic polynomial basis, and $e(b_{i};\gamma)=\textrm{expit}(\gamma_{0}+\gamma_{1}^{\top}b_{i})$. In the importance sampling (IS) estimator proposed in Section 6.1, and in the importance sampling/doubly robust estimator (IS/DR) of Section 6.2, the fitted treatment assignment model was the same, with the outcome model specified through $m(z_{i},s_{i};\phi)=\phi_{0}+\phi_{1}z_{i}+\phi_{2}^{\top}s_{i}$.

For comparison to the Bayesian estimators, we also calculate naive unadjusted comparison (naive), outcome regression adjusted for covariates $s_{i}$ (adjusted), the standard inverse probability of treatment weighted estimator (6) (IPTW), the semi-parametric doubly robust estimator (4) (DR), the clever covariate version of this (CC), the inverse probability of treatment weighted outcome regression based estimator (7) (OR/IPTW), as well as propensity score adjusted outcome regression based estimator (3) (OR/PS). For IPTW, DR, CC, and OR/IPTW, the standard errors were estimated through the standard frequentist nonparametric bootstrap (Efron, , 1979). For OR/PS, to demonstrate the variance estimation issues discussed in Section 4, we calculated both observed information based standard errors, and the adjusted sandwich type standard errors discussed in Supplementary Appendix 3.

The results over 1000 replications are shown in Table 1. Under scenario (I), all estimators except naive and adjusted can correct for confounding, although the joint estimation approach produces a slight bias. In terms of efficiency, the estimators based on propensity score adjusted outcome regression are the best, with the inverse probability of treatment weighting based estimators losing slightly. As discussed in Section 3.2, the clever covariate approach results in higher variability compared to the other doubly adjusted estimators. In terms of variance estimation, the comparison between the unadjusted and adjusted standard errors for OR/PS suggests that under this simulation setting estimation of the propensity scores substantially reduces the variance, and not adjusting for this results in overcoverage. The resampling based variance estimators adjust for this automatically. However, the two-step approach to variance estimation performs poorly; as demonstrated in Supplementary Appendix 3, the two-step point estimator has the same asymptotic variance as the other OR/PS estimators, but the two-step variance estimators unnecessarily add a further variance component. Thus, the simulation results support the discussion in Sections 4 and 5; the two-step and joint estimation approaches are difficult to justify from Bayesian arguments, and do not seem to provide practical advantages in terms of their frequency-based properties. On the other hand, the importance sampling based Bayesian estimators produce very similar results to the OR/IPTW and DR estimators, respectively.

Under scenario (II), all the estimators except IPTW are unbiased, which is expected based on their previously discussed theoretical properties. When the outcome model is correctly specified, there is also very little difference in the efficiency of the various estimators.

8. Discussion

In this paper we reviewed previously proposed Bayesian approaches for propensity score adjusted inferences, focusing on the assumptions concerning correct model specifications. Here it is important to make a distinction between misspecification due to omission of relevant covariates from the outcome model, and misspecification of the functional form of the dependency between the covariates and the outcome. The frequentist propensity score adjusted outcome regression is robust against the former type of model misspecification, but this property is lost in Bayesian estimation, if the misspecified outcome model is allowed to feed back to the estimation of the propensity scores. While feedback issue has been well documented in the literature (e.g. McCandless et al., 2009b, ; Zigler et al., , 2013), and the reasons behind this were already stated by Robins and Ritov, (1997), here we attempted to make the assumptions underlying the Bayesian propensity score approach more explicit. On the other hand, we point out that cutting this feedback in a two-step Bayesian estimation procedure unnecessarily inflates the posterior variance estimates.

As reaching double robustness through Bayesian propensity score adjustment looks difficult, herein we attempted a completely different approach through posterior predictive inferences. Our proposed approach decouples the outcome regression and treatment assignment model through introducing the inverse probability of treatment weights as importance sampling weights in Monte Carlo integration in evaluating posterior predictive expectations. A similar procedure was used in a marginal structural modelling context by Saarela et al., 2015b , improved to its present form in Saarela et al., 2015c . While they used the importance sampling approach for estimating marginal outcome models in a longitudinal setting, herein we showed that in a point treatment setting the combination of importance sampling and posterior predictive inferences can be used to motivate weighted outcome regression or semi-parametric doubly robust inferences. Such a possibility was mentioned, but not formally justified, by Saarela et al., 2015a who applied the importance sampling procedure in the context of estimating optimal treatment regimes. The disadvantage of the importance sampling approach is the same as in the corresponding frequentist inverse probability of treatment weighted inference procedures: the importance sampling weights add variability to the point estimator. In order to control this, a standard approach would be to truncate the weights (e.g. Xiao et al., , 2013), which would also be possible in the importance sampling context (Ionides, , 2008). Recently, Vehtari & Gelman (2015, arXiv:1507.02646v2) suggested probabilistic truncation of importance sampling weights; studying this possibility in the present context is a topic for further research.

Acknowledgement

The authors acknowledge support of the Natural Sciences and Engineering Research Council (NSERC) of Canada.

Supplementary material

Supplementary material available includes Supplementary Appendices 1-4 referred to herein, containing proofs to theorems and other technical material.

Appendix 1

Causal inference as a prediction problem

The estimator (2) can be motivated without the use of potential outcomes notation as a posterior predictive expectation for a new observation under a hypothetical completely randomized setting $\mathcal{E}$ where $Z_{i}{\perp\!\!\!\!\perp_{\mathcal{E}}}X_{i}$ and the probabilities $\mathrm{pr}_{\mathcal{E}}(Z_{i}=a)$ are known constants (cf. the randomized trial measure discussed by Røysland, , 2011). The data are observed under a setting $\mathcal{O}$, where $Z_{i}{{\not\perp\!\!\!\!\perp}_{\mathcal{O}}}X_{i}$, and causal inference then corresponds to inference across these regimes. We can now write for $i\notin\{1,\ldots,n\}$

[TABLE]

where the last form was obtained by choosing the non-parametric specification ${\int_{\psi}p(x_{i};\psi)\pi_{n}(\psi){\,\textrm{d}}\psi}={p(x_{i}\,|\,x_{1},\ldots,x_{n})}={\sum_{k=1}^{n}\delta_{x_{k}}(x_{i})/n}$. Alternatively, in (15) one could use the Bayesian bootstrap (Rubin, , 1981) specification $p(x_{i}\,|\,x_{1},\ldots,x_{n};\xi)=\sum_{k=1}^{n}\xi_{k}\delta_{x_{k}}(x_{i})$, where $\xi=(\xi_{1},\ldots,\xi_{n})$, with $\pi_{n}(\xi)$ taken to be a uniform Dirichlet distribution (see Section 6). Obtaining (16) also required assuming that $p_{\mathcal{E}}(y_{i}\,|\,z_{i},x_{i};\phi)=p_{\mathcal{O}}(y_{i}\,|\,z_{i},x_{i};\phi)\equiv p(y_{i}\,|\,z_{i},x_{i};\phi)$ and and $p_{\mathcal{E}}(x_{i};\psi)=p_{\mathcal{O}}(x_{i};\psi)\equiv p(x_{i};\psi)$, which corresponds to the stability assumption discussed by Dawid and Didelez, (2010).

Appendix 2

Proofs to theorems

Proof** (to Theorem 1).**

If (i) holds true, then also (\mathbf{Y}_{0i},\mathbf{Y}_{1i})\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}e(B_{i})\,|\,S_{i}, and the propensity score adjustment does not add information. If on the other hand (ii) holds true, $\{e(B_{i}),S_{i}\}$ has jointly the balancing property Z_{i}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X_{i}\,|\,\{e(B_{i}),S_{i}\}. This follows from Theorem 2 of Rosenbaum and Rubin, (1983, p. 44) and also implies that (\mathbf{Y}_{0i},\mathbf{Y}_{1i})\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z_{i}\,|\,\{e(B_{i}),S_{i}\} (Rosenbaum and Rubin, , 1983, Theorem 3). Now

[TABLE]

and further,

[TABLE]

*The consistency of the estimator (3) then relies on being able to consistently estimate the quantities in (17). *

Proof** (to Theorem 2).**

Consider first the expectation of (8) under the assumption that (i) holds true. Now

[TABLE]

which followed because now $p(\mathbf{y}_{ai}\,|\,s_{i};\phi)=f(\mathbf{y}_{ai}\,|\,x_{i})$ at the true parameter value. Thus, the misspecified weights do not influence the estimation (in terms of bias) as long as the outcome model is correctly specified.

Under the assumption that (ii) holds true, we have in turn that

[TABLE]

since now $p(z_{i}\,|\,b_{i})=f(z_{i}\,|\,x_{i})$. Using the partitioning $x_{i}=(s_{i},r_{i})$, we can write the last form in above as

[TABLE]

*Thus, even though the outcome regression does not include a sufficient set of confounders, through the IPT weighting we can still obtain valid estimates for the conditional distributions $p(\mathbf{y}_{ai}\,|\,s_{i};\phi)$. *

Proof** (to theorem 3).**

Now the marginal posterior distribution of the parameters $\gamma$ becomes

[TABLE]

Proof** (to Theorem 4).**

The estimator obtained through substituting in the parametric models is

[TABLE]

First, if the outcome model is correctly specified in the sense that $m\{z_{k},s_{k};\widehat{\phi}(\xi)\}=m(z_{k},x_{k};\xi)$, we get

[TABLE]

because the second to last form is equivalent to (The doubly robust estimator as a posterior predictive expectation) in Supplementary Appendix 4. Thus, the first summation term in (Proof) cancels out, leaving only the model based estimator, which itself is now equivalent to the posterior predictive mean difference.

*On the other hand, if the treatment assignment model is correctly specified in the sense that

$\mathrm{pr}_{\mathcal{O}}\{Z_{k}=a\,|\,b_{k};\widehat{\gamma}(\xi)\}=\mathrm{pr}_{\mathcal{O}}(Z_{k}=a\,|\,x_{k},v;\xi)$, we get*

[TABLE]

Therefore, the estimator (Proof) now reduces to

[TABLE]

*which is again equivalent to the posterior predictive mean difference (see Supplementary Appendix 4). *

Appendix 3

On the frequency-based properties of the two-step approach

Trivially, if the outcome model is correctly specified, then \Phi\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}\Gamma\,|\,V and (9) reduces to (16). The interesting situations are naturally those where this is not the case. We denote the log-likelihood by $q_{i}(\phi;\gamma)=\log p\{y_{i}\,|\,z_{i},e(b_{i};\gamma),s_{i};\phi)\}$ and $q(\phi;\gamma)\equiv\sum_{i=1}^{n}q_{i}(\phi;\gamma)$ and consider the quasi-maximum likelihood estimator $\widehat{\phi}(\widehat{\gamma})\equiv\arg\max_{\phi}q(\phi;\widehat{\gamma})$. If the treatment assignment model is correctly specified, $\widehat{\gamma}\rightarrow\gamma_{0}$. In addition, we assume that with any fixed value of $\gamma$, $\widehat{\phi}(\gamma)\rightarrow\phi_{0}(\gamma)$, where $\phi_{0}(\gamma)$ is the parameter vector which minimizes the Kullback-Leibler divergence to the true outcome model (e.g. White, , 1982, p. 4). Thus, by the law of large numbers and continuous mapping, we can write in the usual way that

[TABLE]

where the right hand side is maximized at zero when $\phi=\phi_{0}(\gamma_{0})$, at which point

[TABLE]

Since we also have that the posterior $p(\gamma\,|\,x,y)\rightarrow\delta_{\gamma_{0}}(\gamma)$ in distribution, we can then conjecture that posterior predictive inferences based on (9) will be asymptotically uncounfounded.

With the definitions

[TABLE]

and noting that $U^{\phi}(\widehat{\phi};\gamma^{(j)})=0$ for each $\gamma^{(j)}$, $j=1,\ldots,m$, we can consider the first order Taylor expansion of $U^{\phi}(\widehat{\phi};\gamma^{(j)})$ around the true parameter values $(\phi_{0},\gamma_{0})$, which becomes

[TABLE]

and further,

[TABLE]

if $\frac{1}{m}\sum_{j=1}^{m}\gamma^{(j)}\approx\widehat{\gamma}$. Hence,

[TABLE]

Here we have, by another first order expansion around $\gamma_{0}$, that

[TABLE]

so finally,

[TABLE]

We may similarly expand $U^{\phi}(\widehat{\phi};\widehat{\gamma})$ where the parameters $\gamma$ have been fixed to their maximum likelihood estimates around $(\phi_{0},\gamma_{0})$ as

[TABLE]

to find that

[TABLE]

where

[TABLE]

Since the two estimators $\frac{1}{m}\sum_{j=1}^{m}\widehat{\phi}(\gamma^{(j)})$ and $\widehat{\phi}(\widehat{\gamma})$ have the same linear approximation which is a sum of independent terms, we conclude that they have the same asymptotic distribution,

[TABLE]

Fitting the regression model $y_{i}=\phi_{0}+\phi_{1}z_{i}+\phi_{2}^{\top}s_{i}+\phi_{3}^{\top}g\{e(b_{i};\gamma)\}+\varepsilon_{1i}$ to estimate the parameter of interest $\phi_{1}$ is numerically equivalent to fitting the sequence of regressions $y_{i}=\nu^{\top}s_{i}^{*}(\widehat{\gamma})+\varepsilon_{2i}$, $z_{i}=\alpha^{\top}s_{i}^{*}(\widehat{\gamma})+\varepsilon_{3i}$ and $\{y_{i}-\widehat{\nu}^{\top}s_{i}^{*}(\widehat{\gamma})\}=\phi_{1}\{z_{i}-\widehat{\alpha}^{\top}s_{i}^{*}(\widehat{\gamma})\}+\varepsilon_{4i}$, where $s_{i}^{*}(\gamma)\equiv[s_{i},g\{e(b_{i};\gamma)\}]$. Denoting the estimating function corresponding to the last regression as

[TABLE]

and the partial derivatives of this as e.g. $U^{\phi_{1}\gamma}\equiv\partial U^{\phi_{1}}/\partial\gamma$, we can expand this around $(\phi_{10},\gamma_{0},\nu_{0},\alpha_{0})$, where $\nu_{0}\equiv\nu_{0}(\gamma_{0})$ and $\alpha_{0}\equiv\alpha_{0}(\gamma_{0})$ are the limiting values of the nuisance parameters, as

[TABLE]

since here $\mathrm{E}\{U_{i}^{\phi_{1}\nu}(\phi_{10},\gamma_{0},\nu_{0},\alpha_{0})\}=\mathrm{E}\{U_{i}^{\phi_{1}\alpha}(\phi_{10},\gamma_{0},\nu_{0},\alpha_{0})\}=0$. We can now see that the Theorem 1 of Henmi and Eguchi, (2004, p. 935) applies to the last form here, implying that $\mathrm{avar}(\widehat{\phi}_{1})\leq\mathrm{avar}(\tilde{\phi}_{1})$, where $\widehat{\phi}_{1}$ is the solution to $U(\phi_{1},\widehat{\gamma},\nu_{0},\alpha_{0})=0$ and $\tilde{\phi}_{1}$ is the solution to $U(\phi_{1},\gamma_{0},\nu_{0},\alpha_{0})=0$.

Appendix 4

The doubly robust estimator as a posterior predictive expectation

We first note that because

[TABLE]

we have that

[TABLE]

The usual IPT-weighted estimator for a marginal causal contrast may be derived through a posterior predictive argument as follows. First,

[TABLE]

and thus,

[TABLE]

and

[TABLE]

On the other hand, the usual outcome model based estimator may be motivated similarly as in Supplementary Appendix 1 through

[TABLE]

and

[TABLE]

Finally, we note that we can write (The doubly robust estimator as a posterior predictive expectation) alternatively as

[TABLE]

and therefore

[TABLE]

Thus, the posterior predictive mean difference can be written as

[TABLE]