Bayesian data fusion for unmeasured confounding
Leah Comment Brent A. Coull Corwin Zigler Linda Valeri





















Anharmonicity can enhance the performance of quantum refrigerators
    A. S. Majumdar
    December 30, 2023
==================================================================




Bayesian causal inference offers a principled approach to policy evaluation of proposed interventions on mediators or time-varying exposures. We outline a general approach to the estimation of causal quantities for settings with time-varying confounding, such as exposure-induced mediator-outcome confounders. We further extend this approach to propose two Bayesian data fusion (BDF) methods for unmeasured confounding. Using informative priors on quantities relating to the confounding bias parameters, our methods incorporate data from an external source where the confounder is measured in order to make inferences about causal estimands in the main study population. We present results from a simulation study comparing our data fusion methods to two common frequentist correction methods for unmeasured confounding bias in the mediation setting. We also demonstrate our method with an investigation of the role of stage at cancer diagnosis in contributing to Black-White colorectal cancer survival disparities.






Causal inference; Data fusion; g-formula; Mediation; Racial disparities; Unmeasured confounding.






§ INTRODUCTION

The value of causal evidence from a statistical analysis depends on the quality and suitability of the data source. In the era of Big Data, many large data sources valuable for health research are used in ways that were not foreseen by the original collectors. As such, these sources are missing one or more key covariates. For example, electronic health records and tumor registries may lack important socioeconomic and behavioral factors. If these unmeasured factors act as confounders of the relationship(s) of interest, causal quantities may not be estimated using the observed data, regardless of the sample size. 

Widespread availability of large samples has also ushered in more sophisticated statistical models for disentangling causal effects. Even when many covariates are measured, analyses that take advantage of the rich, longitudinal nature of data sources like electronic health records are vulnerable to time-varying confounding by unmeasured variables. Proper control of confounding is particularly difficult when exposure status changes over time, with later exposure determined in part by covariates influenced by previous exposure. This phenomenon appears in almost every medical setting as doctors tailor treatment based on patient history and current health state. Analogous problems arise in the context of mediation analysis, within which the temporal ordering of treatments, mediators, and outcomes can yield structures analogous to time-varying exposures.

When important confounders are unavailable, researchers typically conduct sensitivity analyses to assess whether bias due to the unmeasured confounding is likely to alter the substantive conclusions of the research. Recent methodological advances have identified sharp nonparametric bounds for common causal estimands such as the average treatment effect <cit.> as well as various mediation quantities <cit.>. Several bias correction formulae provide adjusted point estimates and confidence intervals based on bias values found in the literature <cit.>. In the absence of information about the sources of confounding, they can be used to determine the strength of confounding needed to eliminate statistical significance. Individual approaches also require rare outcomes, specific link functions, or assumptions about effect modification <cit.>. With some notable exceptions <cit.>, correction methods rarely incorporate uncertainty surrounding the bias parameters. Generally speaking, most existing sensitivity methods suffer from poor extensibility to both arbitrary confounding structures and longitudinal settings.

Fortunately, the era of Big Data is also the era of abundant data. Relationships among the outcome, confounders, and exposure of interest can be found in alternative data sets, though these sources may not be as representative of the target population as the main source. A literature on data fusion methods has arisen to meet the need to combine information from multiple sources. Recent authors have proposed Bayesian variable selection methods with validation data sets <cit.> and data integration for information from different scales, such as individual-level and community-level data  <cit.>.

To address the limited extensibility and uncertainty quantification of existing methods, we propose a general framework for incorporating information from external data sources to perform sensitivity analyses for unmeasured confounding. We develop a Bayesian method for data fusion and incorporate it into an existing Bayesian g-formula approach which adjusts for confounding using parametric models for covariate standardization <cit.>. To handle unmeasured confounding, we introduce a procedure for generating informative priors using external data sources. We then describe two estimation strategies to account for unmeasured confounding: one using mechanics similar to <cit.>, and the other augmenting this approach with a Bayesian bootstrap procedure for marginalization. We compare such strategies with traditional sensitivity analysis approaches, paying particular attention to potential violations of causal transportability <cit.> when the underlying causal processes differ between the two populations from which the main and external data were sampled. Unlike existing approaches, these Bayesian g-formula methods generalize to accommodate unmeasured confounding of many different types, including time-varying confounding found in mediation and analysis of longitudinal treatments.

 Our paper is organized as follows. Section 2 develops a Bayesian g-method that accommodates dynamic and stochastic treatment assignment mechanisms, then highlights the connection between this model and a mediation analysis. We describe two estimators based on the g-formula. Section 3 introduces two Bayesian data fusion algorithms to implement sensitivity analyses for unmeasured confounding. A simulation study comparing the two methods to traditional sensitivity analysis approaches is given in Section 4. In Section 5, we use the data fusion method from Section 3 to augment cancer registry data with information from a cohort study in order to evaluate the role of stage at diagnosis in explaining Black-White disparities in colorectal cancer survival. We conclude with a discussion in Section 6.
 


§ THE BAYESIAN G-FORMULA FOR STATIC AND DYNAMIC REGIMES WITHOUT UNMEASURED CONFOUNDING



 §.§ Causal notation and assumptions

Let Y denote the observed outcome of interest in a causal graph G, with the central scientific question concerning two intervention regimes g and g'. One of these regimes may correspond to the “natural” assignment mechanism that generated the observed data. Let V be the intervention set (i.e., the nodes intervened upon by either g or g'), and let Z be the set of baseline confounders and post-treatment variables not influenced by treatment. Let W be the set containing all other nodes in G, in which case W includes any variable influenced by treatment but not of primary interest (i.e., not the outcome Y) or directly intervened upon (i.e., W ∉V). The complete set of observed data is O = (Z, V, W, Y). Let Y^g denote the potential outcome for Y under regime g, with the causal contrast of interest τ = Y^g - Y^g'. Depending on the specifics of the regime, mediators of the A → Y relationship may either be in the set W or V. Potential outcomes for W and V under regime g are denoted by W^g and V^g. When one mediator is of primary interest, we will denote it by M. Throughout this paper, we must assume Bayesian analogs to positivity (Assumption 1) and consistency (Assumption 2) <cit.>. We also require exchangeability conditional on observed variables (Assumption 3) and correct specification of all parametric models (Assumption 4). Formal statements of these assumptions can be found in the supplemental materials.



 §.§ Static and deterministic treatment regimes

After adopting parametric models indexed by the parameter vector θ, the Bayesian g-formula algorithm outlined by <cit.> gives the posterior predictive distribution for a newly observed outcome Y under intervention regime g_0 ∈{g, g' } as

    p(ỹ^g_0 | o ) = ∫ p(ỹ^g_0 | θ, o) π(θ|o) dθ

where π(θ|o) is the posterior of the parameters θ given the observed data O. 
The posterior distribution of the causal effect τ is therefore

    p(τ| o) = ∫( p(ỹ^g | θ, o) - p(ỹ^g' | θ, o) ) π(θ|o) dθ .

Keil and colleagues outline a simulation-based algorithm for estimating causal contrasts for static regimes. To facilitate our extension to the mediation setting, we introduce different notation to emphasize the distinction between covariates Z that are unaffected by treatment and covariates W that are caused by one or more variables in the intervention set. For parametric models p(z | θ_Z), p(w| v, z, θ_W), and p(y | w, v, z, θ_Y), the respective parameter likelihoods are ℒ(θ_Z | z), ℒ(θ_W | w, v, z), and ℒ(θ_Y | y, w, v, z). Then the likelihood for the complete parameter vector θ = (θ_Z, θ_W, θ_Y) is given by

    ℒ(θ|o) = 
    ℒ(θ_Y | y, w, v, z)
    ×ℒ(θ_W | w, v, z) 
    ×ℒ(θ_Z | z).

We assume that θ_Y, θ_W, θ_Z are independent a priori such that π(θ) = π(θ_Y) ×π(θ_W) ×π(θ_Z). The resulting parameter posterior is π(θ|o) ∝ℒ(θ|o) ×π(θ), and the posterior predictive distribution for the outcome under regime g_0 ∈{g, g'} is

    p(ỹ^g_0 | o) =    
    
    
    p(ỹ | g_0, w̃, z̃, θ_Y)
    p(w̃ | g_0, z̃, θ_W)
    p(z̃ | θ_Z)
    
    
    π(θ|o)
    
    dθ
    dw̃
    dz̃.


To distinguish Keil's computationally intensive approach from classical covariate standardization techniques using the g-formula, we refer to this approach as simulation-based Bayesian g-formula (BGF-SIM). We now develop a procedure to accommodate dynamic and stochastic treatment regimes.



 §.§ Dynamic and stochastic treatment regimes

Scientific questions of interest sometimes involve contrasts of regimes which assign exposure stochastically according to different distributions depending on prior covariates. In particular, we may be interested in the “natural” assignment mechanism generating the observed data. This is exactly the case for mediation analysis, which decomposes the effect of an exposure on an outcome into component causal pathways in order to understand possible mechanisms enacting the overall effect. Figure <ref> shows a classic causal structure in mediation, where the mediator M channels part of the effect of the exposure A on the outcome Y, with Z as a baseline confounder. The set W contains other mediators which also act as exposure-induced mediator-outcome confounders for the M → Y relationship. 



All common mediation estimands can be formulated as contrasts in regimes, including controlled direct effects (CDE), natural direct and indirect effects (NDE and NIE), and randomized interventional analogs to these quantities <cit.>. For concreteness, we restrict attention in the main text to the randomized interventional analog to the natural direct effect (rNDE), which is identified under weaker conditions than the NDE; estimation algorithms for the CDE, NDE, NIE, and rNIE are available in the supplemental materials. The rNDE can be conceived as a contrast in dynamic regimes where part of the regime recreates the naturally occurring assignment mechanism <cit.>. Because the naturally occurring assignment mechanism for M is unknown, it must be estimated. We add a parametric model for M reflecting its parents in the natural stochastic regime, and let θ = (θ_Z, θ_W, θ_M, θ_Y). The likelihood conditional on observed data is

    ℒ(θ|o) = ℒ(θ_Y | y, m, w, a, z)
    ×ℒ(θ_M|m, w, a, z)
    ×ℒ(θ_W | w, a, z)
    ×ℒ(θ_Z | z).

We now outline a strategy for estimating the rNDE which compares a=1 to a=0 when the mediator is stochastically assigned as it would be under a=0. The rNDE contrasts regime g=(A=1, M=H_z(a=0)) with g'=(A=0, M=H_z(a=0)), where H_z(a=0) is a draw from the distribution of M under A=0, conditional on Z. That is, p(h_z(a=0)|z̃) = ∫ p(m̃|a=0, z̃, w̃)p(w̃|a=0,z̃)dw̃.

The equation for p(ỹ^g' | o) for g'=(A=0, M=H_z(a=0)) is analogous to Equation <ref> with an added model for M. For the regime g=(A=1, M=H_z(a=0)), the W value used to assign M is different from the value for Y (i.e., a recanting witness) <cit.>. The simulation-based Bayesian g-formula resolves this issue with an independence assumption resulting in separate posterior predictive draws of w̃^a for both a=0 and a=1. The posterior mean of the rNDE is thus given by

    rNDE = [    ỹ( 
    p(ỹ | a=1, m, w̃^1, z̃, θ_Y)
    p(w̃^1 | a=1, z̃, θ_W) - 
    
       
    p(ỹ | a=0, m, w̃^0, z̃, θ_Y)
    p(w̃^0 | a=0, z̃, θ_W)) ×
       
    p(m̃ | a=0, z̃, w̃^0, θ_M) 
    p(z̃ | θ_Z)
    π(θ|o)
    ]
    dθ
    dw̃^0
    dw̃^1
    dm̃
    dz̃
    dỹ.




 §.§ The Bayesian g-formula with the Bayesian bootstrap: an alternative marginalization strategy for closed-form estimands

Previous sections assumed that marginalization over the baseline confounder distribution occurred through posterior predictive sampling from p(z̃|θ_Z). In practice, Z can be high-dimensional, and parametrically modeling p(z|θ_Z) introduces additional opportunities for model misspecification. Because Z is by definition the same for all regimes, <cit.> suggest sampling z̃ nonparametrically from the observed empirical distribution of Z, p_N(z). For settings where Z, A, and M are all discrete, we introduce an alternative approach. Because the g-formula estimators have closed forms in the discrete case, marginalization can occur through the Bayesian bootstrap. The Bayesian bootstrap assigns observation weights (d_1, …, d_n) sampled from a Dirichlet(1,…,1), with weights changing every MCMC iteration. The closed-form causal contrast τ is calculated for every observed Z, yielding τ(Z_i) for i=1,…,n. Then the weighted average ∑_i=1^n d_i τ(Z_i) gives the posterior draw of the population average causal effect for that MCMC iteration. We refer to this procedure as the closed-form Bayesian g-formula (BGF-CF).

For relatively simple causal graphs with discrete data, BGF-CF avoids the computationally intensive posterior predictive simulation of BGF-SIM. Many popular causal estimands have tractable closed form solutions for discrete data. As an example, the rNDE for baseline confounder value z, rNDE(z), can be written as 

    ∑_w, m{   𝔼[Y|a=1,w,m,z] p(w|a=1,z) 
    - 
    𝔼[Y|a=0,w,m,z] p(w|a=0,z)
    }
    p(m|a=0,z).

Model-based estimates of every term in Equation <ref> can obtained with each posterior parameter sample θ^(b), with p(m|a=0,z) = ∑_w' p(m|a=0,w',z,θ_M^(b))p(w'|a=0,z,θ_W^(b)). If there are K distinct baseline confounder patterns and B MCMC iterations, the Bayesian g-formula only requires calculating rNDE(z) a total of B × K times. Letting ξ_k denote the sample frequency of covariate pattern k, the Bayesian bootstrap weights (d_1,…,d_K) are repeatedly sampled from a Dirichlet(ξ_1,…,ξ_K). At each iteration, the posterior draw of the population rNDE is ∑_k=1^K d_k × rNDE(z_k)^(b). For K ≪ n, this method can be much more computationally efficient than the simulation-based g-formula. 




§ BAYESIAN DATA FUSION FOR UNMEASURED CONFOUNDING

We now consider the problem of making causal inferences when an important confounder is unmeasured in the primary (“main”) data set but is available in a secondary (“external”) source. Although the data fusion algorithm we outline accommodates arbitrary confounding structures and many different estimands, we restrict attention to estimating a randomized natural direct effect with an exposure-induced unmeasured mediator-outcome confounder U as in Figure <ref>. This setting is interesting for two reasons: (1) the additional complexity involved with estimating the natural, stochastic assignment mechanism and (2) the need to accommodate exposure-induced confounding. Complete estimation algorithms for the rNIE, CDE, and average treatment effects for longitudinal exposures are all available in the supplemental materials.





 §.§ Notation and assumptions

Suppose that an investigator is interested in learning about an effect in some target population for which a large data source (N=n_1) exists. The desired causal quantity is the rNDE in the population from which the n_1 observations were randomly sampled. Figure <ref> shows the causal structure, with outcome Y, exposure A, mediator M, and baseline confounders Z. There is also an exposure-induced mediator-outcome confounder U, which is unmeasured in the data set of size n_1. Information on {A,M,Z,Y,U} exists in a smaller secondary data source (N = n_2, with n_2 < n_1), which may or may not be from the same population. For g-computation, models for {A, M, Y} can be any univariate or multivariate generalized linear model, and there are no distributional restrictions on Z. However, imposing additional restrictions on U can dramatically improve MCMC performance. Marginalization over the distribution of U is only guaranteed for U with finite support, although certain continuous distributions can also be integrated out of the likelihood.

In addition to the assumptions of the Bayesian g-formula <cit.> – which are formally stated as Assumptions 1-4 in the Web Appendix – Bayesian data fusion requires causal transportability to hold. That is, although the smaller data set may not be representative of the target population with respect to the distribution of baseline covariates, the underlying causal processes operate in the same way.

(Parametric causal transportability)
Let 𝒫_1 and 𝒫_2 denote the superpopulations for the main and external data sources. If 𝒫_1 𝒫_2, the causal graph structures of 𝒫_1 and 𝒫_2 must agree such that all X ∈{U, M, Y} have the same parent nodes pa(X). Furthermore, the true underlying data-generative parameters θ_X, 𝒫_1 and θ_X, 𝒫_2 must be the same such that

    p(x | pa(X), θ_X, 𝒫_1) = 
    p(x | pa(X), θ_X, 𝒫_2).

Note that causal transportability holds by design if the external data are a random validation sample from 𝒫_1.


The conditional exchangeability of Assumption 3 is still required for U (i.e., U^a  A | Z, a conditional independence which is encoded in Figure <ref>), but the requirements for M and Y are relaxed to allow the confounder U to be unmeasured. 
(Conditional exchangeability)
Briefly, (Z,U) must be sufficient to control confounding. For the randomized rNDE in Figure <ref>, this implies:

    M^a  A     | Z, U 
    
    Y^a  A     | Z, U, M





 §.§ Specification of parametric models


To illustrate the closed-form estimator and facilitate contrasts with existing methods, we assume A, U, M, and Y are all binary with logistic link functions. The baseline confounders Z are also assumed to be discrete. 

Letting π_X = P(X = 1| pa(X)), we adopt the following models: 

    = 
    γ_0 + γ_A A_i + z_i'γ_Z
     = 
    β_0 + β_A A_i + z_i'β_Z + β_U U_i 
     =  
    α_0 + α_A A_i + z_i'α_Z + α_M M_i + α_AMA_i M_i + α_U U_i

Then θ_U = (γ_0, γ_A, γ_Z), θ_M = (β_0, …, β_U), and θ_Y = (α_0, …, α_U). For X ∈{U, M, Y}, 
(θ_X|x, pa(X)) = f(x|pa(X), θ_X) = ∏_i=1^N (π_X,i)^x_i (1-π_X,i)^1-x_i. Equation <ref> shows the observed data likelihood for the full parameter vector θ = (θ_U, θ_M, θ_Y) in the main data set, marginalizing over the missing U.

    ℒ_m = ∏_i=1^n_1[ ∑_uℒ(θ_Y | y_i, m_i, u_i=u, a_i, z_i)
    ℒ(θ_M|m_i, u_i=u, a_i, z_i)
    ℒ(θ_U|u_i=u, a_i, z_i)
    ]


For a generic prior π(θ), the posterior for θ marginalizing over the missing U is proportional to ℒ_m ×π(θ). How to set an informative prior π(θ) using the secondary data set is the focus of the next section.



 §.§ Specification of prior information with external data
 
Given that U is unmeasured in the main data source, any parameters involving U (i.e., θ_U, β_U, and α_U) cannot be identified from that data. Because the main data set is presumably more representative of the target population of interest, the sole reason for integrating the external data set is for providing information about the confounder unmeasured in the main data set. That information can be summarized through the use of informative priors.

Priors for {θ_U, θ_M, θ_Y } are derived by fitting in the external data frequentist maximum likelihood models corresponding to Equations <ref> through <ref>. Under causal transportability, maximum likelihood estimators fit in the external data will be consistent and asymptotically normal about θ_𝒫_1. For X∈{U,M,Y}, let θ̂_X,MLE denote the maximum likelihood estimate (MLE) of θ_X in the external data, and let Σ_X,MLE be the estimated variance-covariance matrix. Then 𝒩(θ̂_X,MLE, Σ_X,MLE) is a sensible choice for π(θ_X). With moderately large n_2, this prior approximates the posterior distribution for θ in a Bayesian analysis conducted using the n_2 observations, assuming a non-informative prior. With a priori independence, the complete prior for θ is  π(θ)=π(θ_U)×π(θ_M)×π(θ_Y).

If 𝒫_1 𝒫_2, then less informative priors may be preferable for the identifiable parameters. Consider the parameter α_A, about which the main data source contains substantial information. In the model formulation given by Equation <ref>, variance and covariance hyperparameters for α_A would be found along the second row and column of Σ̂_Y, MLE. If we multiply the off-diagonal elements in that row and column by a large inflation factor (e.g., σ=1000) and the diagonal element by σ^2, we assert a marginal prior distribution on α_A that is virtually non-informative. However – critically –  the prior correlation between α_A and the unidentifiable parameter α_U is preserved.



 §.§ A simulation-based Bayesian data fusion algorithm (BDF-SIM)

We now outline a simulation-based Bayesian data fusion approach for g-formula causal contrasts in the context of rNDE estimation.


  * Fit maximum likelihood models in the external data to obtain the prior π(θ) as detailed in Section <ref>.

  * Use a No-U-Turn sampler (NUTS) with target probability distribution proportional to ℒ_m ×π(θ) in order to obtain posterior samples of the regression parameter vector θ. The probabilistic programming language Stan has a NUTS implementation <cit.>, and it is available to R users through the  R package <cit.>. For some large B (e.g., 4,000), let  θ^(1), …, θ^(B) denote the B posterior samples remaining after discarding warmup iterations.

  * For MCMC iteration b=1,…,B and i = 1,…, n_1:


  * Sample baseline covariate vector z̃_i from the empirical distribution.

  * For g_0 ∈{ g, g' } and a_0 ∈{ 0 , 1 }, sample ũ^a_0, g_0 (b)_i as Bernoulli with success probability 

    γ_0^(b) + z̃_i' γ_Z^(b) + a_0 γ_A^(b)


  * For g_0 ∈{ g, g' }, sample randomized mediator m̃^0, g_0 (b)_i as Bernoulli with success probability

    β_0^(b) + z̃_i' β_Z^(b) + β_U^(b)ũ^0, g_0 (b)_i


  * Define individual-level causal contrast ϕ̃_i^(b) as 

    ϕ̃_i^(b) =    α_0^(b) + z̃_i' α_Z^(b) + α_M^(b)m̃^0, g (b)_i + α_A^(b) + α_AM^(b)m̃^0, g (b)_i +
    α_U^(b)ũ^1, g (b)_i
        -  
    α_0^(b) + z̃_i' α_Z^(b) + α_M^(b)m̃^0, g' (b)_i +
    α_U^(b)ũ^0, g' (b)_i



  * Calculate population estimate rNDE^(b) = ∑_i=1^n_1ϕ̃_i^(b) / n_1.

  * Construct a point estimate for rNDE as the posterior mean rNDE=∑_b=1^B rNDE^(b)/B, and create quantile-based 95% credible intervals as the 2.5^th and 97.5^th quantiles of (rNDE^(1), …, rNDE^(B)).




 §.§ A Bayesian data fusion algorithm for closed-form estimands using the Bayesian bootstrap (BDF-CF)

We now outline a data fusion procedure for the rNDE using the Bayesian g-formula with the Bayesian bootstrap for confounder marginalization. Steps 1, 2, and 5 of the closed-form version with the Bayesian bootstrap are identical to the simulation-based approach, so we show only steps 3 and 4.



  * For b=1,…,B and k=1,…,K for the K unique covariate patterns


  * Sample covariate patterns weights (d_1^(b), …, d_K^(b)) from a Dirichlet(ξ_1,…,ξ_K), where ξ_k is the count of observations with the unique covariate pattern Z_k.

  * Calculate the individual-level contrast ϕ(z_k) for pattern z_k according to Equation <ref>, replacing W with U and plugging in the appropriate model-based estimates from θ^(b). Concretely,

    𝔼[Y|a,u,m,z] 
    =    α_0^(b) + z' α_Z^(b) + α_M^(b) m + α_A^(b) a +
    α_AM^(b) am +
    α_U^(b) u 
    
    
    p(u|a,z) 
    =    γ_0^(b) + γ_A^(b) a + z' γ_Z^(b)
    
    
    p(m|a=0,z)
    =    ∑_u(
    β_0^(b) + z' β_Z^(b) + β_U^(b) u p(u|a=0,z) )



  * Calculate population estimate rNDE^(b) = ∑_k=1^K (d_k^(b)×ϕ_k^(b))/n_1.


Although the simulation-based and closed-form Bayesian g-formula approaches are identical with respect to regression parameter estimation, their differences have implications for extensibility to other causal estimands and scalability to large data sets. The performance of these two estimators under various conditions is the focus of the next section.




§ SIMULATION STUDY


We designed a simulation study to evaluate the performance of the simulation-based and closed-form BDF estimators relative to existing bias corrections we briefly describe in this section. The estimand of interest was the rNDE in the main study superpopulation. Due to the fact that sensitivity analyses based on sharp nonparametric bounding and those based on externally derived bias parameters are not directly comparable, we focus our comparison between BDF and other bias correction techniques.
 


 §.§ Data generation procedure

We considered a number of scenarios with varying data generation schemes. We varied: sample sizes (n_1 and n_2 =n_1/10), causal structure (Δ_U,A = 1 if the mediator-outcome confounder U is caused by A and 0 otherwise), and presence of an interaction (Δ_Y,AM = 1 if there is an A-M interaction in the Y model and 0 otherwise). The strength of mediator-outcome confounding by U was governed by two quantities, β_U and α_U, the log-odds ratios of U in the M and Y models, respectively. When the same β_U and α_U were used to generate the main and external data sets, we have causal transportability; this was done for strong confounding by U (β_U = α_U = 1.5, corresponding to odds ratios of ≈ 4.5). To investigate the performances of the various approaches under violations of the transportability assumption, β_U = α_U = 0 was used to generate the external data, while β_U = α_U = 1.5 in the main data. Complete details of the data generation process are available in the Web Appendix. For each simulation condition, estimator bias and interval coverage were assessed using 200 replicates.



 §.§ Implementation of the Bayesian data fusion estimators

Closed-form and simulation-based variants of Bayesian data fusion were implemented for each pair of simulated main and external data sets. Priors were constructed using the external data as described in Section <ref> without variance inflation for the covariance matrix (i.e., assuming transportability). The models in the external data were correctly specified, with Δ_Y,AM and Δ_U,A matching the underlying generation process for the main data. Posterior samples of the bias-corrected rNDE were obtained from 3 MCMC chains of 2,000 iterations each, with the first 1,000 samples discarded as warmup. The posterior mean was taken as a point estimate, with uncertainty captured using 95% quantile-based credible intervals.



 §.§ Alternative bias correction methods

The first comparator method, referred to as the delta-gamma (DG) correction, is a classical bias correction method <cit.>. A version for controlled direct effects can be used for the rNDE when the two coincide, i.e., if (1) U is not exposure-induced, and (2) there is no exposure-mediator interaction in the outcome model. Note that for estimands on the risk difference scale, (2) does not hold for logistic models of Y even if the A-M interaction coefficient is zero. This approach also requires that the effect of U should be the same across all levels of A (i.e., Y|a,z,m,U=1 - Y|a,z,m,U=0 does not depend on a), which cannot hold in a logistic model unless the A coefficient is zero. For comparability with BDF we elected to use the secondary data source, replacing component quantities in the DG bias formula with maximum likelihood estimates derived from logistic regression models in the external data. Confidence intervals were obtained by bootstrapping the main data 200 times. 

A second frequentist correction, referred to as the interaction correction (IX), can accommodate exposure-mediator interaction in the outcome model <cit.>. Originally derived as a bias correction for the NDE, it is more generally applicable than the DG correction, but it similarly requires that U not be exposure-induced. Again, we fit maximum likelihood models to the external data source to derive bias-corrected estimates within each covariate pattern, and confidence intervals were obtained via the bootstrap. Additional details on the DG and IX  implementation are available in the supplemental materials.



 §.§ Simulation results

Figure <ref> shows estimates from the case where U is exposure-induced. When the transportability assumption holds, both BDF estimators eliminate the confounding bias at all sample sizes. In contrast, the frequentist correction methods do worse than no correction at all. Although these corrections do not purport to address exposure-induced mediator-outcome confounding, this finding underscores the danger of using these corrections when U may be caused by A.

In the absence of transportability, the information extracted from the external data set by the BDF procedure is misleading, and the estimators perform poorly. Confounding bias is not eliminated, and the prior information leads to less posterior uncertainty. The frequentist estimators also do not correct the bias, but the uncertainty is the same as the uncorrected naive intervals.



Table <ref> reports coverage for the 95% confidence and credible intervals; the Web Appendix contains information regarding interval widths. In general, credible intervals from BDF approaches had widths comparable to the naive and frequentist corrected confidence intervals. However, the classical correction methods are not unbiased, and the interval coverage is low (< 10%). Conversely, as noted previously, BDF methods did not perform well in the case of transportability violations, achieving less than 1% coverage.





§ EXAMINING THE ROLE OF STAGE AT DIAGNOSIS IN BLACK-WHITE SURVIVAL DISPARITIES IN COLORECTAL CANCER



 §.§ Overview

We now use BDF to explore the extent to which differentials in stage at diagnosis contribute to apparent racial disparities in colorectal cancer survival. Our analysis provides an estimate of how much we could reasonable expect to reduce the observed survival disparity if we could break the  between race and cancer stage at the time of diagnosis, e.g., by implementing targeted screening programs that lead to earlier colorectal cancer detection among Blacks.

<cit.> sought to address this question in a recent article with data from a registry of US cancer patients from 1992-2005. The National Cancer Institute's Surveillance, Epidemiology, and End Results (SEER) registry collects information on tumor site and stage for a sizable proportion of cancer patients from diverse geographic regions within the US. They concluded that eliminating Black-White disparities in colorectal cancer stage at diagnosis would lead to a 35% reduction in survival disparities as measured by 5-year restricted mean survival time. Their analysis controlled for a number of covariates, including age at diagnosis, gender, time period of cancer diagnosis, geographic locale, and median county income as derived from the American Community Survey; however, it did not control for household-level poverty status, as that information was not available.



 §.§ Analysis description

We extend the analysis of Valeri and colleagues by incorporating information about confounding of the stage-survival relationship by individual-level income using data from the Cancer Care Outcomes Research and Surveillance (CanCORS) Consortium data. This observational study followed patients shortly after cancer diagnosis and aimed “to determine how the characteristics and beliefs of lung and colorectal cancer patients, physicians and health-care organizations influence treatments and outcomes spanning the continuum of cancer care from diagnosis to recovery or death, and to evaluate the effects of specific therapies on patients’ survival, quality of life, and satisfaction with care” <cit.>. As a result of these ambitious aims, the CanCORS database contains detailed socioeconomic information, including household income for the year preceding cancer diagnosis. We chose U=1 to correspond the lowest income group of <$40,000 per year. The goal was to assess the bias of the residual disparity measure as calculated in SEER, assuming true underlying race-poverty and poverty-survival relationships in SEER matched those estimated in CanCORS. The survival outcome was a binary indicator Y for whether the patient was alive 5 years post-cancer diagnosis. Self-reported race was coded such that A=1 for non-Hispanic blacks and A=0 for non-Hispanic whites; individuals reporting Hispanic origin were excluded. The intervening variable of interest, stage at cancer diagnosis M, took on values 1-4 corresponding to cancer stages I-IV. Adjustment covariates included in all models were: gender, age at cancer diagnosis (<60, 60-65, or >65), and geographic region (West, South, or other). Patients whose cancer was unstaged were excluded, leaving a total of 146,031 colorectal cancer cases in the SEER analysis data set.

First, we fit two naive models using maximum likelihood in the SEER data: (1) stage at cancer diagnosis as a function of race and adjustment covariates, using a baseline category logit model; and (2) 5-year survival as a function of race, stage at diagnosis, and the adjustment covariates, using a logistic link. Coefficients from these models were used to calculate a naive estimated residual disparity measure RD_naive and bootstrapped 95% confidence intervals.

Next, we implemented both BDF estimators to obtain poverty-adjusted estimates of the black-white survival disparity. To construct priors, we fit three frequentist models using the 1,613 CanCORS colorectal cancer patients for whom complete stage and covariate data were available. The two regression models described above were modified by adding poverty as a covariate. Since SEER is more representative of the target population for intervention, all parameters except the bias coefficients were given marginally noninformative prior distributions with the variance inflation strategy outlined in Section <ref>. A third and final frequentist model was a logistic regression for poverty as a function of race, gender, region, and age category. Because none of these U-related parameters are identifiable in SEER, no variance inflation procedure was performed.

Using simulation and closed-form BDF, we estimated the poverty-adjusted residual disparity in the SEER data. Four chains of 2,000 MCMC iterations each were run in Stan <cit.>, with the first 1,000 iterations discarded as warmup. The Gelman-Rubin convergence diagnostic R̂ was calculated for all parameters <cit.>.



 §.§ Residual disparity results

Posterior samples of the poverty-adjusted population residual disparity measure calculated using simulation-based and closed-form BDF are shown in Figure <ref>. The null value of zero, which represents Black-White equality with respect to baseline covariate-adjusted survival, lies beyond the far right of the graph. Visible as a dotted line on the left is the disparity we currently observe without an intervention on stage. With a value of 0.099 (95% CI: 0.092, 0.107), we estimate that Black patients are 9.9 percentage points less likely to survive 5 years post-diagnosis than White patients of the same gender and geographic region. The naive estimate of the residual disparity after an intervention aligning Blacks' cancer stage distribution to the current stage suggests that the remaining disparity in 5-year survival would be 6.6% (95% CI: 5.8, 7.4). The BDF analyses suggest that unmeasured confounding by poverty does not substantially change the estimated residual disparity, with closed-form and simulation-based estimates of 6.5% (95% CI: 5.3, 7.4).

Given the abundance of literature documenting the role of socioeconomic status in cancer outcomes <cit.>, it may be surprising to see adjustment for poverty having such a small impact on the estimated residual disparity. One possibility is that causal transportability may not hold between the SEER and CanCORS populations in ways related to poverty (i.e., which cannot be addressed with variance inflation). That is, CanCORS may appear to be representative of the larger US population from which SEER draws its cancer cases <cit.>, but the causal relationships determining cancer outcomes in CanCORS are fundamentally different from the processes in SEER because many CanCORS study sites are academic medical centers in large cities <cit.>. Thus, we may not see dramatic shifts in our conclusions for SEER because CanCORS does not contain evidence for substantial stage-survival confounding by poverty. Alternatively, there may residual confounding due to the coarsening of socioeconomic deprivation – a complex, multifaceted problem – into a single binary indicator. Nevertheless, this analysis integrating the two data sources gives policymakers two potentially valuable pieces of information: (1) a quantitative estimate of the poverty-adjusted residual disparity and (2) a better understanding of the true uncertainty surrounding that estimate. 





§ DISCUSSION

In this paper, we have proposed a general method for Bayesian data fusion that can be used to perform sensitivity analyses for unmeasured confounding in a variety of settings. The approach addresses forms of confounding in both static and dynamic treatment regimes as well as in mediation, including exposure-induced mediator-outcome confounding. While there can be no substitute for a well-designed study in the target population of interest, decision makers cannot wait for the ideal analysis in the ideal data set and must often work from incomplete or imperfect information. Bayesian data fusion communicates the sensitivity of a research conclusion while incorporating some of what is already known about the problem.

The general nature of the BDF-SIM and BDF-CF data fusion methods make them easily extendable to other settings. For example, any number of parametric models could be used for the unmeasured confounder. Although we demonstrated properties using the randomized natural direct effect in a mediation setting, these principles can be applied to any mediation estimand or to settings with time-varying confounding. With BDF-SIM, any generalized linear model can be adopted for the unmeasured confounder, allowing for both continuous and discrete distributions. Both BDF-SIM and BDF-CF can accommodate multiple unmeasured confounders. Depending on the types of confounding present, information on multiple confounders could be constructed from different external data sources, although doing so may require some assumptions about the joint distribution of their effects in the outcome model. 

With respect to the motivating question of Black-White racial disparities in U.S. colorectal cancer patients, we conclude that unmeasured confounding of the stage-survival relationship by poverty leads to residual disparity reduction estimates that are slightly too optimistic. Implementing an intervention – for example, a targeted screening program – to alleviate or eliminate delayed cancer diagnosis for Black colorectal cancer patients would substantially improve 5-year survival outcomes. However, without also intervening upon the complex societal factors that lead to greater poverty among black patients, we cannot realize the full benefit of such an intervention.

Both the closed-form and simulation-based variations of BDF are limited in part by their reliance on two major assumptions. First, they assume that the parametric models are otherwise correctly specified, which means they do not account for uncertainty in model misspecification, as some sensitivity analyses do <cit.>. A more comprehensive uncertainty quantification would incorporate additional uncertainty due to model selection. Second, as with most parametric causal inference, the models extrapolate causal effects <cit.>, and problems with overlap must be detected by the analyst. Recent advances in Bayesian nonparametrics <cit.> may be adapted to add flexibility to portions of the models. Third, like other data fusion methods, BDF assumes causal transportability between the external and main data sets <cit.>. We can perform exploratory analyses investigating representativeness, but we can never be certain that the underlying causal generative processes in two study populations are comparable. With BDF, the skeptical analyst can modify the variance-covariance matrices in the prior distribution to have greater variances, effectively reducing the “prior sample size” of the information on the unmeasured confounder. However, doing so places greater prior probability mass on confounding parameter values that may be implausible. Taken to the extreme, the external data source provides no information at all. When it comes to transportability and data fusion, one must balance caution about non-transportability and the desire to let the secondary data source inform the choice of sensitivity parameters. Thus, for data-driven sensitivity analysis methods to gain widespread acceptance, the statistical literature needs more formal methods for evaluating and characterizing causal transportability.

Nevertheless, these BDF methods represent significant steps forward in statistical analyses to inform decisionmaking. As a sensitivity analysis method, the flexible nature of BDF-SIM allows for a much greater variety of causal structures and regression model specifications, and it works for any causal estimand that can be represented as a counterfactual contrast. As a data fusion approach, the underlying Bayesian principles allow for extensibility to handle uncertainty quantification for multiple unmeasured confounders of various types. Information from the external data source enters exclusively through prior distributions, reducing computational burden and sidestepping data privacy concerns. BDF only needs  parameter estimates and variance-covariance matrices from the external data set; sharing this information does not compromise the privacy of individual study participants. This fact dramatically increases the number of data sources which may be used as external data. Researchers in the era of Big Data cannot guarantee that the right data are always available to them, but developing statistical methods to rigorously synthesize information from multiple sources can give decision-makers the tools to make more informed choices.






§ ACKNOWLEDGEMENTS

LC acknowledges support from NIH grants T32CA009337 and T32ES007142 as well as a David H. Peipers Fellowship. CZ was supported by EPA grant RD-835872 and NIH grants R01ES026217 and NIH R01GM111339. LV received funding from NIH grant K01MH118477, and BC was funded by NIH grants P01CA134294 and P30ES000002.








[Antonelli, Zigler, and DominiciAntonelli
  et al.2017]antonelli2017guided
Antonelli, J., Zigler, C., and Dominici, F. (2017).
Guided bayesian imputation to adjust for confounding when combining
  heterogeneous data sources in comparative effectiveness research.
Biostatistics page kxx003.

[Avin, Shpitser, and PearlAvin
  et al.2005]avin2005identifiability
Avin, C., Shpitser, I., and Pearl, J. (2005).
Identifiability of path-specific effects.
Department of Statistics, UCLA .

[Ayanian, Chrischilles, Wallace, Fletcher,
  Fouad, Kiefe, Harrington, Weeks, Kahn, Malin, et alAyanian
  et al.2004]ayanian2004understanding
Ayanian, J., Chrischilles, E., Wallace, R., Fletcher, R., Fouad, M., Kiefe, C.,
  Harrington, D., Weeks, J., Kahn, K., Malin, J., et al. (2004).
Understanding cancer treatment and outcomes: the cancer care outcomes
  research and surveillance consortium.
Journal of Clinical Oncology 22, 2992–2996.

[Carpenter, Gelman, Hoffman, Lee, Goodrich,
  Betancourt, Brubaker, Guo, Li, and RiddellCarpenter
  et al.2016]carpenter2016stan
Carpenter, B., Gelman, A., Hoffman, M., Lee, D., Goodrich, B., Betancourt, M.,
  Brubaker, M., Guo, J., Li, P., and Riddell, A. (2016).
Stan: A probabilistic programming language.
Journal of Statistical Software 20,.

[Catalano, Ayanian, Weeks, Kahn, Landrum,
  Zaslavsky, Lee, Pendergast, and HarringtonCatalano
  et al.2013]catalano2013representativeness
Catalano, P., Ayanian, J., Weeks, J., Kahn, K., Landrum, M. B., Zaslavsky, A.,
  Lee, J., Pendergast, J., and Harrington, D. (2013).
Representativeness of participants in the cancer care outcomes
  research and surveillance (CanCORS) consortium relative to the surveillance, epidemiology and end results (SEER) program.
Medical Care 51, e9.

[Didelez, Dawid, and GenelettiDidelez
  et al.2012]didelez2012direct
Didelez, V., Dawid, P., and Geneletti, S. (2012).
Direct and indirect effects of sequential treatments.
arXiv preprint arXiv:1206.6840 .

[Ding and VanderWeeleDing and
  VanderWeele2016]ding2016sensitivity
Ding, P. and VanderWeele, T. (2016).
Sensitivity analysis without assumptions.
Epidemiology 27, 368.

[Ding and VanderweeleDing and
  Vanderweele2016]ding2016sharp
Ding, P. and Vanderweele, T. (2016).
Sharp sensitivity bounds for mediation under unmeasured
  mediator-outcome confounding.
Biometrika 103, 483–490.

[Gelman and RubinGelman and
  Rubin1992]gelman1992inference
Gelman, A. and Rubin, D. (1992).
Inference from iterative simulation using multiple sequences.
Statistical science pages 457–472.

[GreenlandGreenland2005]greenland2005multiple
Greenland, S. (2005).
Multiple-bias modelling for analysis of observational data.
JRSS:A 168, 267–306.

[Jackson, Best, and RichardsonJackson
  et al.2006]jackson2006improving
Jackson, C., Best, N., and Richardson, S. (2006).
Improving ecological inference using individual-level data.
Statistics in Medicine 25, 2136–2159.

[Keil, Daza, Engel, Buckley, and EdwardsKeil
  et al.2015]keil2015bayesian
Keil, A., Daza, E., Engel, S., Buckley, J., and Edwards, J. (2015).
A bayesian approach to the g-formula.
Statistical Methods in Medical Research page 0962280217694665.

[Le, Ziogas, Lipkin, and ZellLe
  et al.2008]le2008effects
Le, H., Ziogas, A., Lipkin, S., and Zell, J. (2008).
Effects of socioeconomic status and treatment disparities in
  colorectal cancer survival.
Cancer Epidemiology and Prevention Biomarkers 17,
  1950–1962.

[McCandless and SomersMcCandless and
  Somers2017]mccandless2017bayesian
McCandless, L. and Somers, J. (2017).
Bayesian sensitivity analysis for unmeasured confounding in causal mediation analysis.
Statistical Methods in Medical Research page 0962280217729844.

[Pearl and BareinboimPearl and
  Bareinboim2011]pearl2011transportability
Pearl, J. and Bareinboim, E. (2011).
Transportability of causal and statistical relations: A formal
  approach.
In Data Mining Workshops (ICDMW), 2011 IEEE 11th International
  Conference on, pages 540–547. IEEE.

[Pearl and BareinboimPearl and
  Bareinboim2014]pearl2014external
Pearl, J. and Bareinboim, E. (2014).
External validity: From do-calculus to transportability across
  populations.
Statistical Science pages 579–595.

[Roy, Lum, Zeldow, Dworkin, Re III, and
  DanielsRoy et al.2017]roy2017bayesian
Roy, J., Lum, K., Zeldow, B., Dworkin, J., Re III, V. L., and Daniels, M.
  (2017).
Bayesian nonparametric generative models for causal inference with
  missing at random covariates.
Biometrics .

[Stan Development TeamStan Development
  Team2016]rstan
Stan Development Team (2016).
RStan: the R interface to Stan.
R package version 2.14.1.

[Tchetgen Tchetgen and ShpitserTchetgen
  Tchetgen and Shpitser2012]tchetgen2012semiparametric
Tchetgen Tchetgen, E. and Shpitser, I. (2012).
Semiparametric theory for causal mediation analysis: efficiency
  bounds, multiple robustness, and sensitivity analysis.
Annals of Statistics 40, 1816.

[Valeri, Chen, Garcia-Albeniz, Krieger,
  VanderWeele, and CoullValeri et al.2016]valeri2016role
Valeri, L., Chen, J., Garcia-Albeniz, X., Krieger, N., VanderWeele, T., and
  Coull, B. (2016).
The role of stage at diagnosis in colorectal cancer Black–White
  survival disparities: a counterfactual causal inference approach.
Cancer Epidemiology and Prevention Biomarkers 25, 83–89.

[VanderWeeleVanderWeele2010]vanderweele2010bias
VanderWeele, T. (2010).
Bias formulas for sensitivity analysis for direct and indirect
  effects.
Epidemiology 21, 540.

[VanderWeeleVanderWeele2015]vanderweele2015explanation
VanderWeele, T. (2015).
Explanation in causal inference: methods for mediation and
  interaction.
Oxford University Press.

[VanderWeele and ChibaVanderWeele and
  Chiba2014]vanderweele2014sensitivity
VanderWeele, T. and Chiba, Y. (2014).
Sensitivity analysis for direct and indirect effects in the presence
  of exposure-induced mediator-outcome confounders.
Epidemiology, Biostatistics, and Public Health 11,.

[Vansteelandt, Bekaert, and
  ClaeskensVansteelandt et al.2012]Vansteelandt2012uncertainty
Vansteelandt, S., Bekaert, M., and Claeskens, G. (2012).
On model selection and model misspecification in causal inference.
Statistical Methods in Medical Research 21, 7–30.
PMID: 21075803.






§ SUPPLEMENTAL MATERIALS







§ ASSUMPTIONS FOR THE BAYESIAN G-FORMULA



[Positivity]
Let X be any node in the causal graph G, and x any value in the support of X. Then for any regime g_0 ∈{ g, g'}, p(x| pa_g_0(X) = p̃ã) > 0, where p̃ã is any value taken by pa_g_0(X), the parent nodes of X under g_0. Furthermore, it must hold for all x and p̃ã that p(x| pa(X) = p̃ã) > 0, where pa(X) without a subscript indicates the parent nodes of X in the naturally occurring treatment assignment mechanism.

[Consistency]
For any regime g_0 ∈{g, g'}, Y = Y^g_0 whenever V takes on the values prescribed by g_0. If V is a single binary treatment, this statement simplifies to Y = V Y^1 + (1 - V)Y^0.

[Conditional exchangeability]
For any variable V_0 in the intervention set and every regime g_0 ∈{g, g' }, there exists a set of measured variables C ⊂{Z, W} such that Y^g_0 V_0 | C. 


[Correct parametric model specification]
For every node X ∈{V, W, Y} modeled conditional on variables C with parameters θ_X, the parametric model f(X | C,  θ_X) is correctly specified.

 


§ BAYESIAN G-FORMULA ALGORITHMS FOR OTHER CAUSAL ESTIMANDS

For simplicity of exposition, we assume a single unmeasured confounder U throughout the following algorithms, but U can be vector-valued, or there may be L distinct unmeasured confounders U_1, U_2, …, U_L occupying different positions in the causal graph. Although the algorithms outlined below are the simulation-based BDF-SIM, the more computationally efficient BDF-CF is available when all variables are discrete and in selected other instances (e.g., certain outcome models with identity link functions). Throughout, we continue using n_1 to denote the sample size of the main data sources.



 §.§ Time-varying confounding of a longitudinal exposure



Suppose that the true causal DAG is as in Figure <ref>, with a discrete exposure-induced unmeasured variable U acting as a confounder of the time-varying exposure A measured at two time points to yield A = (A_1, A_2). Denote the regimes of interest with g = (A_1 = a_1, A_2 = a_2) and g' = (A_1 = a_1', A_2 = a_2'), with the causal estimand of interest being the superpopulation average causal effect ACE = Y^g - Y^g'.


Adopt parametric generalized linear models indexed by θ_X for X ∈{U, A_2, Y}, with h_X(·) denoting the link function and η_X the linear predictor term, which is a function only of the parent nodes X and θ_X. Equation <ref> gives a general model form. 

    h_X(X_i | pa(X_i), θ_X) =    η_Xi

The intervention set is V = (A_1, A_2). For g_0 ∈{g, g'} let pa^g_0(X) be the intervened parents of X where any V ∈{pa(X) } has been set (deterministically or stochastically) in accordance with g_0 and the relevant model.

To emphasize that η_Xi depends on the parameters θ_X as well as the values of the parents of X for observation i, we can also write it as η_Xi(pa(X), θ_X). 

Let the  ℒ(θ_X | x , pa(X)).
Then θ = (θ_U, θ_A_2, θ_Y) and

    ℒ_c =    ∏_i=1^n_1ℒ(θ_Y | y_i, a_2i, u_i, a_1i, z_i)
    ℒ(θ_A_2 | a_2i, u_i, a_1i, z_i)
    ℒ(θ_U | u_i, a_1i, z_i)

For discrete U, this yields the marginal likelihood of

    ℒ_m =    ∏_i=1^n_1[ ∑_u
    ℒ(θ_Y | y_i, a_2i, u_i = u, a_1i, z_i)
    ℒ(θ_A_2 | a_2i, u_i = u, a_1i, z_i)
    ℒ(θ_U | u_i = u, a_1i, z_i) ]




  * Fit maximum likelihood models in the external data to obtain the prior π(θ) as detailed in Section <ref>.

  * Use NUTS with target probability distribution proportional to ℒ_m ×π(θ) in order to obtain posterior samples of the regression parameter vector θ. For some large B, let  θ^(1), …, θ^(B) denote the B posterior samples remaining after discarding warmup iterations.

  * For MCMC iteration b=1,…,B and i = 1,…, n_1:


  * Sample baseline covariate vector z̃^(b)_i from the empirical distribution.

  * For g_0 ∈{g, g' }, set ã_̃1̃^g_0 (b)_i deterministically or stochastically in accordance with g_0. For example, if g is the static, deterministic regime setting A_1 to level a_1, ã_̃1̃^g (b)_i = a_1 for all i and b.

  * For g_0 ∈{g, g' }, sample ũ^g_0 (b)_i in accordance with

    h_X^-1( η_X( pa^g_0(x̃^g_0_i), θ_X^(b)) )

Concretely, for g = (A_1 = a_1, A_2 = a_2) and a logistic model P(U_i = 1 | Z_i, A_1i) = γ_0 + γ_A_1 A_1i + γ_Z'Z_i, sampling ũ^g (b)_i requires drawing from a Bernoulli with success probability

    γ_0^(b) + γ_A_1 a_1 + γ_Z^(b)'z̃^(b)_i



  * For g_0 ∈{g, g' }, set ã_̃2̃^g_0 (b)_i deterministically or stochastically in accordance with g_0. Concretely, if g is the static, deterministic regime setting A_2 to level a_2, ã_̃2̃^g (b)_i = a_2 for all i and b.

  * For g_0 ∈{g, g' }, draw ỹ^g_0 (b)_i in accordance with g_0, θ_Y^(b), z̃^(b)_i, ũ^g_0 (b)_i, ã_̃1̃^g_0 (b)_i, and ã_̃2̃^g_0 (b)_i. Then calculate the individual-level causal effect

    ϕ̃_i^(b) =    ỹ^g (b)_i - ỹ^g' (b)_i

Alternatively, if the conditional mean of Y has a closed form μ(θ_Y, z, u, a_1, a_2), define the individual-level causal effect as

    ϕ̃_i^(b) =    μ(θ_Y^(b), z̃^(b)_i, ũ^g (b)_i, ã_̃1̃^g (b)_i, ã_̃2̃^g (b)_i ) - 
    μ(θ_Y^(b), z̃^(b)_i, ũ^g' (b)_i, ã_̃1̃^g' (b)_i, ã_̃2̃^g' (b)_i ).

For example, if Y_i is conditionally normal with mean α_0 + α_Z' Z_i + α_A_1 A_1i + α_U U_i + α_A_2 A_2i + α_intx× A_1i× A_2i× U_i and the contrast of interest compares regimes g = (A_1 = a_1, A_2 = a_2) and g' = (A_1 = a_1', A_2 = a_2'), 

    ϕ̃_i^(b) = 
       α_A_1^(b) a_1 + α_U^(b)ũ^g (b)_i + α_A_2^(b) a_2 + α_intx^(b) a_1 × a_2 ×ũ^g (b)_i 
     - 
       α_A_1^(b) a_1' + α_U^(b)ũ^g' (b)_i + α_A_2^(b) a_2' + α_intx^(b) a_1' × a_2' ×ũ^g' (b)_i



  * Calculate population estimate ACE^(b) = ∑_i=1^n_1ϕ̃_i^(b) / n_1.

  * Construct a point estimate for ACE as the posterior mean ACE=∑_b=1^B ACE^(b)/B, and create quantile-based 95% credible intervals as the 2.5^th and 97.5^th quantiles of (ACE^(1), …, ACE^(B)).




 §.§ Natural direct effects

Suppose the true causal diagram underlying the mediation is shown in Figure <ref>. The unmeasured confounder U can confound the (1) exposure-mediator, (2), exposure-outcome, or (3) mediator-outcome relationships, as well as any combination of (1) - (3).


The population average natural direct effect of changing exposure A to a instead of a^*, while holding the mediator M to its natural value under A = a^*, is given by NDE = Y^a M^a^* - Y^a^* M^a^*. This estimand has an intervention set V = { A, M } and can be formulated as a contrast in the regimes g = (A = a, M = M^a^*) and g' = (A = a^*, M = M^a^*).



  * Fit maximum likelihood models in the external data to obtain the prior π(θ) as detailed in Section <ref>.

  * Use NUTS with target probability distribution proportional to ℒ_m ×π(θ) in order to obtain posterior samples of the regression parameter vector θ. For some large B, let  θ^(1), …, θ^(B) denote the B posterior samples remaining after discarding warmup iterations.

  *  For MCMC iteration b=1,…,B and i = 1,…, n_1:


  * Sample baseline covariate vector z̃^(b)_i from the empirical distribution.

  * Sample ũ^(b)_i in accordance with θ_U^(b) and z̃_i^(b) using

    h_U^-1( η_U( z̃_i^(b), θ_U^(b)))

In contrast to previous algorithms, this sampling does not need to be done for each g_0 ∈{ g, g' } because U cannot be a descendant of A or M in the causal graph for the natural direct effect to be well defined.

  *  Sample mediator m̃^a^* (b)_i according to θ_M^(b), z̃_i^(b), and ũ^(b)_i using

    h_M^-1( η_M( z̃_i^(b), ũ^ (b)_i, a^*, θ_M^(b)))



  * For g_0 ∈{g, g' }, draw ỹ^g_0 (b)_i in accordance with g_0, θ_Y^(b), z̃^(b)_i, ũ^(b)_i, and m̃^a^* (b)_i. For regime g = (A = a, M = M^a^*) this involves 

    h_Y^-1( η_Y( z̃_i^(b), ũ^ (b)_i, m̃^a^* (b)_i, a, θ_Y^(b)))

while for g' = (A = a^*, M = M^a^*) the relevant equation will involve 

    h_Y^-1( η_Y( z̃_i^(b), ũ^ (b)_i, m̃^a^* (b)_i, a^*, θ_Y^(b)))

Then calculate the individual-level causal effect

    ϕ̃_i^(b) =    ỹ^g (b)_i - ỹ^g' (b)_i

Alternatively, if the conditional mean of Y has a closed form μ(θ_Y, z, u, a, m), define the individual-level causal effect as

    ϕ̃_i^(b) =    μ(θ_Y^(b), z̃^(b)_i, ũ^(b)_i, a, m̃^(b)_i ) - 
    μ(θ_Y^(b), z̃^(b)_i, ũ^(b)_i, a^*, m̃^a^* (b)_i ).



  * Calculate population estimate NDE^(b) = ∑_i=1^n_1ϕ̃_i^(b) / n_1.

  * Construct a point estimate for NDE as the posterior mean NDE=∑_b=1^B NDE^(b)/B, and create quantile-based 95% credible intervals as the 2.5^th and 97.5^th quantiles of (NDE^(1), …, NDE^(B)).




 §.§ Natural indirect effects

Again suppose that the correct causal diagram is as in Figure <ref>, where natural direct effects are well defined.

The population average natural indirect effect is the effect of changing the mediator M from the value it naturally takes under exposure A = a^* to the value it naturally takes under A=a, while holding the exposure constant at level a. In potential outcome notation, this quantity is given by NIE = Y^a M^a - Y^a M^a^*. This estimand has an intervention set V = { A, M } and can be formulated as a contrast in the regimes g = (A = a, M = M^a) and g' = (A = a, M = M^a^*).

The estimation algorithm is the same as in Section <ref> until Step <ref>, where it continues as follows.

enumi3

  * For MCMC iteration b=1,…,B and i = 1,…, n_1:


  * Sample baseline covariate vector z̃^(b)_i from the empirical distribution.

  * Sample ũ^(b)_i in accordance with θ_U^(b) and z̃_i^(b) using

    h_U^-1( η_U( z̃_i^(b), θ_U^(b)))


  * For each a_0 ∈{a, a^*}, sample mediator m̃^a_0, (b)_i according to θ_M^(b), z̃_i^(b), and ũ^(b)_i using

    h_M^-1( η_M( z̃_i^(b), ũ^ (b)_i, a_0, θ_M^(b)))



  * For g_0 ∈{g, g' }, draw ỹ^g_0 (b)_i in accordance with g_0, θ_Y^(b), z̃^(b)_i, ũ^(b)_i, and m̃^a_0 (b)_i for the a_0 corresponding to the M counterfactual in g_0. For regime g = (A = a, M = M^a) this involves 

    h_Y^-1( η_Y( z̃_i^(b), ũ^ (b)_i, m̃^a (b)_i, a, θ_Y^(b)))

while for g' = (A = a, M = M^a^*) the relevant equation will involve 

    h_Y^-1( η_Y( z̃_i^(b), ũ^ (b)_i, m̃^a^* (b)_i, a, θ_Y^(b)))

Then calculate the individual-level causal effect

    ϕ̃_i^(b) =    ỹ^g (b)_i - ỹ^g' (b)_i

Alternatively, if the conditional mean of Y has a closed form μ(θ_Y, z, u, a, m), define the individual-level causal effect as

    ϕ̃_i^(b) =    μ(θ_Y^(b), z̃^(b)_i, ũ^(b)_i, a, m̃^a (b)_i ) - 
    μ(θ_Y^(b), z̃^(b)_i, ũ^(b)_i, a, m̃^a^* (b)_i ).



  * Calculate population estimate NIE^(b) = ∑_i=1^n_1ϕ̃_i^(b) / n_1.

  * Construct a point estimate for NIE as the posterior mean NIE=∑_b=1^B NIE^(b)/B, and create quantile-based 95% credible intervals as the 2.5^th and 97.5^th quantiles of (NIE^(1), …, NIE^(B)).




 §.§ Randomized interventional analogs to the natural indirect effect

The estimation algorithm is the same as in Section <ref> until Step <ref>, where it continues as follows.
enumi3

  * For MCMC iteration b=1,…,B and i = 1,…, n_1:


  * Sample baseline covariate vector z̃_i from the empirical distribution.

  * For each g_0 ∈{ g, g' } and a_0 ∈{ a, a^*}, sample ũ^a_0, g_0 (b)_i in accordance with θ_U^(b) and z̃_i^(b) using

    h_U^-1( η_U( z̃_i^(b), a_0, θ_U^(b)))



  * For g_0 ∈{ g, g' }, sample randomized mediator m̃^g_0 (b)_i in accordance with θ_U^(b), z̃_i^(b), and ũ_i^a_0, g_0(b). For regime g = (A = a, M = H_z(a = a)), draw m̃^g (b)_i using

    h_M^-1( η_M( z̃_i^(b), a, ũ_i^a, g (b),  θ_M^(b)))

and for g' = (A = a, M = H_z(a = a^*)) draw m̃^g' (b)_i using

    h_M^-1( η_M( z̃_i^(b), a^*, ũ_i^a^*, g' (b),  θ_M^(b)))


  * Define individual-level causal contrast For g_0 ∈{g, g' }, draw ỹ^g_0 (b)_i in accordance with g_0, θ_Y^(b), z̃^(b)_i, ũ^(b)_i, and m̃^g_0 (b)_i. For regime g = (A = a, M = H_z(a = a)), draw ỹ^g (b)_i using

    h_Y^-1( η_Y( z̃_i^(b), ũ^a, g (b)_i, m̃^g (b)_i, a, θ_Y^(b)))

while for g' = (A = a, M = H_z(a = a^*)), draw ỹ^g' (b)_i using

    h_Y^-1( η_Y( z̃_i^(b), ũ^a, g' (b)_i, m̃^g' (b)_i, a, θ_Y^(b)))

Then calculate the individual-level causal effect

    ϕ̃_i^(b) =    ỹ^g (b)_i - ỹ^g' (b)_i

Alternatively, if the conditional mean of Y has a closed form μ(θ_Y, z, u, a, m), define the individual-level causal effect as

    ϕ̃_i^(b) =    μ(θ_Y^(b), z̃^(b)_i, ũ^a, g (b)_i, a, m̃^g (b)_i ) - 
    μ(θ_Y^(b), z̃^(b)_i, ũ^a, g' (b)_i, a, m̃^g' (b)_i ).



  * Calculate population estimate rNIE^(b) = ∑_i=1^n_1ϕ̃_i^(b) / n_1.

  * Construct a point estimate for rNIE as the posterior mean rNIE=∑_b=1^B rNIE^(b)/B, and create quantile-based 95% credible intervals as the 2.5^th and 97.5^th quantiles of (rNIE^(1), …, rNIE^(B)).




 §.§ Controlled direct effects

The controlled direct effect is the effect of changing exposure A to level a from a^* while holding the mediator M fixed at level m, i.e., CDE = Y^am - Y^a^*m. In general, this requires no unmeasured exposure–outcome confounding and  no unmeasured mediator–outcome confounding. Both of these cases can be addressed by with slight modifications to previously stated versions of the BDF-SIM algorithm in order to obtain B posterior samples CDE^(1), …, CDE^(B).



  §.§.§ Exposure-outcome confounding or mediator-outcome confounding that not affected by treatment

Suppose the true causal diagram is as in Figure <ref>, where U acts as an exposure-outcome confounder, mediator-outcome confounder, or both. (It may also be an exposure-mediator confounder, but the controlled direct effect is already identified if both the U → A and U → Y arrows are missing.)

For this causal structure, the controlled direct effect can be estimated using the algorithm from Section <ref>, replacing the stochastic assignment of m̃^a^* (b)_i in Step <ref> with universal assignment to m for all i and b. 



  §.§.§ Exposure-induced mediator-outcome confounding

Suppose the true causal diagram is as in Figure <ref> in the main text, i.e., where U is an exposure-induced mediator-outcome confounder. Further suppose scientific interest lies in the controlled direct effect of changing A to level a from a^* while holding the mediator M fixed at level m, i.e., CDE = Y^am - Y^a^*m.

The controlled direct effect can be estimated using the algorithm from Section <ref>, replacing A_1 with A and A_2 with M. The two regimes are g = (A = a, M = m) and g' = (A = a^*, M = m).



§ DATA GENERATION PROCEDURE FOR SIMULATIONS

The “no-interaction” case corresponding to no statistical interaction has Δ_Y,AM = 0.

When U is not exposure-induced, Δ_U,A = 0.

For violations of transportability, β_U = α_U = 0 was used for generation of the small data set. Otherwise, β_U = α_U = 1.5 in order to induce strong mediator-outcome confounding by U.

    Z_1    ∼Bernoulli(0.5) 
    
    Z_2 | Z_1    ∼Bernoulli(0.5) 
    
    A | Z_1, Z_2    ∼Bernoulli( -0.2 + 0.5 Z_1 + 0.7 Z_2)
    
    U | A, Z_1, Z_2    ∼Bernoulli( -0.4 + Δ_U,A 1.5 A)
    
    M | U, A, Z_1, Z_2    ∼Bernoulli( -1.5 + 0.3 Z_1 + 0.2 Z_2 + 0.7 A + β_U U)
    
    Y | M, U, A, Z_1, Z_2    ∼Bernoulli( -2 + 0.3 Z_1 + 0.2 Z_2 + A + 0.8 M + Δ_Y, AM A M + α_U U)




§ FREQUENTIST BIAS CORRECTIONS



 §.§ δ-γ correction

For each level of z, the bias due to U for the rNDE on the difference scale, assuming A does not cause U is

    B^CDE_dg,add(m=0|z) =    δ_m=0(z) γ_m=0(z)  with
    δ_m=0(z) =     P(U = 1 | z, m=0, a = 1) - P(U = 1 | z, m=0, a = 0) 
    γ_m=0(z) =    Y=1| z, a, m=0, u = 1 - Y=1| z, a, m=0, u = 1

The DG-corrected population estimate of the rNDE is then given by

    rNDE_dg = ∑_z( rNDE_uc(z) - B^CDE_dg,add(m=0|z) ) p(z).





 §.§ Interaction correction

The bias in the additive NDE from by a mediator-outcome confounder U which is not exposure-induced is given by:

    B^NDE_ix,add(z) 
    =    ∑_m,u( 
    [ Y|a,m,z,uP(u|a,m,z) - Y|a^*,m,z,uP(u|a^*,m,z) ] P(m|a^*, c)
    ) 
        - ∑_m,u(
    [ Y|a,m,z,u - Y|a^*,m,z,u]P(m|a^*,z,u)P(u|z)
    )

In our context, a = 1 and a^* = 0. To eliminate the possibility of model misspecification in the estimation of the bias correction term B^NDE_ix,add(z), saturated parametric logistic regression models were adopted for U and M when possible. (For the data application, sparseness in the covariates made this impossible, and parametric models with many interaction terms were fit to reduce, but not eliminate, misspecification.) The model for Y used to obtain Y|a=1,m,z,u was the same logistic model used to obtain the naive rNDE, except with U as an additional term.

The corrected estimate of the population rNDE was calculated by

    rNDE_ix = ∑_z( rNDE_uc(z) - B^NDE_ix,add(z) ) p(z).




§ ADDITIONAL SIMULATION RESULTS



 §.§ Credible interval widths


Mean confidence and credible interval widths for the simulated scenarios are given in Table <ref>. The Bayesian estimators tend to have widths comparable to frequentist analogs when there is no transportability, but tends to have wider intervals when there is substantial bias in the external data. BDF-CF may perform slightly better than BDF-SIM (i.e., have narrower intervals) in small samples due to smaller Monte Carlo error.





 §.§ Performance when unmeasured confounder is not exposure-induced


In the absence of transportability, all of the bias corrections considered perform poorly. However, the BDF-SIM and BDF-CF estimates are more confident about their (incorrect). In the top right panel of Figure <ref>, one can see that the BDF approaches reduce bias better than the frequentist corrections in smaller sample sizes, but that the difference in performance is virtually eliminated with a main data set size of n=10,000. Interestingly, the IX correction does not exhibit good performance when there is a strong exposure-mediator interaction, particularly in small samples.






