Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
78,866 characters · 24 sections · 68 citation commands
Bayesian Model Averaging in Causal Instrumental Variable Models
{\it Keywords:} Conditional Bayes factors, Endogeneity, Invalid instruments, Identification, Unobserved confounding
Instrumental variables (IV) offer a way to infer causal effects in the presence of unobserved confounding. To be suitable, an instrumental variable must meet two criteria: It must not be affected by the unobserved confounder and must not affect the outcome directly (valid), and it must be associated with the regressor of interest (relevant). In practice, finding variables that fulfill these assumptions is challenging. Furthermore, when multiple instrumental variables are available, the results can be highly sensitive to the specific instruments chosen. In many applications, there is a large degree of uncertainty over which instruments to use in addition to the usual uncertainty over covariate inclusion in a regression model. To address this uncertainty, we propose a Bayesian model averaging (BMA) procedure that averages across different sets of instrumental variables and covariates in a multivariate sampling model with correlated residuals. In our context, endogenous variables are treatments that are affected by unobserved confounding.
An important contribution of this paper is a data-driven approach for selecting valid and relevant instruments. The user specifies a matrix of potential instruments and covariates, and the BMA framework identifies valid and relevant instruments as those variables included in the treatment model but excluded from the outcome model. This approach provides robustness against the use of invalid or irrelevant instruments. We establish a consistency result showing that the posterior asymptotically concentrates on the correct model.
We build on a large existing literature. Key ideas and methodology for Bayesian variable selection in multivariate regression models can be found in Brown_etal_98, while rossi_bayesian_2005 propose a simple Bayesian approach to instrumental variables based on conjugate normal and inverse Wishart priors and a Gibbs sampler iterating between the outcome, treatment, and covariance parameters. koop_bayesian_2012 were among the first to propose BMA in IV models, offering a general approach for simultaneous equations. However, their complex reversible jump MCMC scheme (which can mix quite slowly in practice) makes their method less accessible for practitioners. karl_instrumental_2012 build on the Bayesian IV framework of rossi_bayesian_2005, using conditional Bayes factors to iteratively update the outcome and treatment models within a Gibbs sampler. Subsequent work lenkoski_sovereign_2020,kourtellos2020measuring,lee_incorporating_2022 extends their approach to allow for non-Gaussian and multiple endogenous variables. lenkoski_two-stage_2014 introduce a hybrid Bayesian version of the classical Two-Stage Least Squares (TSLS) estimator.
Our proposed methodology follows karl_instrumental_2012, and we take close inspiration from their computational strategy: we use a Gibbs sampler to update the outcome, treatment, and covariance parameters separately and use conditional Bayes factors in the model updates. However, we also provide a careful analysis of prior structures, and we use suitably adapted $g$-priors on the regression coefficients, which make the analysis scale-invariant, and a flexible Cholesky-based prior on the structural covariance lopes_bayesian_2014. Our approach uses independent priors on the outcome and treatment model space without enforcing identification through zero prior probability on non-identified models. We accommodate non-Gaussian outcomes and endogenous variables by adapting the latent Gaussian framework with univariate link functions (ULLGM) proposed in steel_model_2024. This greatly extends the applicability of our method while not substantially adding to the computational cost. We prove that our model selection procedure is consistent in the sense that the conditional Bayes factors used in the model updates tend to infinity in favor of the true model as the sample size grows. Finally, we provide software in the accompanying gIVBMA.jl package written in the Julia language. While karl_instrumental_2012 previously released an R package for their method, it has since been withdrawn from CRAN. We refer to the method of karl_instrumental_2012 as IVBMA and call our proposed procedure generalised IVBMA (gIVBMA).
Related ideas and methods can also be found in the context of Bayesian inference on binary treatment effects under endogenous treatment selection. Methodology for these models with Bayesian variable selection was developed by jacobi_etal_16 using a single latent factor, and by wagner_etal_23 from a bi-factor model perspective.
An important strand of the classical literature focuses on identification and estimation when some instrumental variables are invalid ditraglia_using_2016, kang_instrumental_2016, windmeijer_use_2019, windmeijer_confidence_2021. kang_instrumental_2016 show that classical identification is still possible as long as a plurality rule holds, {i.e. }the valid instruments outnumber the invalid ones. Our approach performs well in settings with potentially invalid instruments, but is not bound by the plurality rule.
Our proposed method also performs well in cases with many individually weak instruments. Traditional classical estimators tend to exhibit substantial bias and can even become inconsistent in such settings. To address these issues, common approaches include first-stage regularisation\footnote{Inspired by the TSLS procedure, “first-stage” refers to the estimation of the treatment model, while “second-stage” relates to inference in the outcome model.} okui_instrumental_2011, belloni_sparse_2012, carrasco_regularization_2012, using jackknife-fitted values in the first stage angrist_jackknife_1999, hansen_instrumental_2014, selecting instruments based on minimizing a mean-square error criterion donald_choosing_2001, or model-averaging to obtain model-averaged first-stage predictions of the endogenous variable kuersteiner_constructing_2010.
Through extensive experiments on real and simulated data, we evaluate gIVBMA against IVBMA, BMA, and a number of classical methods tailored for this problem. The gIVBMA methodology works well in settings with potentially invalid instruments and many individually weak instruments. Another advantage of gIVBMA (and IVBMA) over most classical methods is that we obtain the posterior distribution of the structural covariance matrix which informs us on the degree of endogeneity. This also provides these methods with the ability to borrow strength from the model for the endogenous variables.
Section {(ref)} introduces the sampling model and discusses model uncertainty and identification. Sections {(ref)} and (ref) respectively describe the prior specification and posterior inference, while Section (ref) provides the model consistency results. Simulated and real data are analysed in Sections {(ref)} and (ref). Finally, Section {(ref)} concludes.
We consider a structural model for $n$ observations of a single outcome and $l\ge 1$ potentially endogenous variables:
where $\vec{y}$ is an $n \times 1$ outcome vector, $\vec{X} = \vec{[x_1 : \ldots : x_l]}$ is an $n \times l$ matrix of causes or treatments that are potentially endogenous, $\vec{Z}$ is an $n \times p$ matrix of instruments and exogenous covariates (i.e. $\vec{Z}\perp\!\!\!\!\perp \vec{\epsilon}$), $\vec{\epsilon}$ is an $n \times 1$ vector of outcome residuals, and $\vec{H = [\eta_1 : \ldots : \eta_l]}$ is an $n \times l$ matrix of treatment residuals. The outcome model is parameterized by an intercept $\alpha$, an $l \times 1$ vector of “effects” $\vec{\tau}$, and a $p \times 1$ vector of instrument (or covariate) coefficients $\vec{\beta}$. The treatment model is parameterized by a $1 \times l$ (row) vector $\vec{\Gamma}$ and a $p \times l$ matrix of slope coefficients $\vec{\Delta}$. The main parameter of interest is the vector of treatment effects $\vec{\tau}$.
We allow all instruments to enter the outcome model, which is essential for our model selection procedure to identify valid and relevant instruments in a data-driven manner. For valid instruments, $\vec{\beta}$ will have most of its posterior mass on (or very close to) zero, but we do not enforce this a priori. This point will become clearer in the next section, where we introduce model uncertainty.
This setup can be motivated by the potential outcomes framework imbens_instrumental_2014. Suppose that $y_i({x})$ is the potential outcome for the $i$-th outcome observation with conditional expectation $ \mathbb{E} \left[ y_i({x}) \mid \vec{Z_i} \right] = \alpha + {x} {\tau} + \vec{Z_i \beta} $ such that the average treatment effect (ATE) of increasing $x$ by one unit is exactly $\tau$. If we have multiple endogenous variables, i.e. $\vec{x}$ is a vector, any component of $\vec{\tau}$ is the ATE of only increasing the corresponding component of $\vec{x}$ by 1 while the other components remain constant. Define the outcome residual $\epsilon_i = y_i(x) - \mathbb{E} \left[ y_i(x) \mid \vec{Z_i} \right]$, which is by definition uncorrelated with the exogenous variables $\vec{Z_i}$. Given observations $(y_i, \vec{X}_i, \vec{Z}_i)$, this potential outcomes model implies the observed data model $ y_i = y_i(\vec{X}_i) = \alpha + \vec{X}_i \vec{\tau} + \vec{Z}_i \vec{\beta} + \epsilon_i, $ which is our outcome equation for a single observation. Thus, $\vec{\tau}$ has a causal interpretation. Subsection (ref) contains more discussion on the identification of $\vec{\tau}$.
We assume the residuals follow a matrix normal distribution (defined in supplementary Section (ref)) $ \vec{[\epsilon : H]} \sim MN(\vec{0}, \vec{I_n}, \vec{\Sigma})$, where the structural covariance matrix $\vec{\Sigma}$ can be partitioned into
If $\vec{\Sigma_{yx}} \neq 0$, the outcome and the treatment residuals are correlated, signaling the presence of unobserved confounding or endogeneity. We assume the errors to be homoskedastic and serially uncorrelated across observations but this could be generalised relatively easily.
The conditional distribution of $\vec{y} \mid \vec{X}$ is given by (see supplementary Section (ref))
where $\sigma_{y|x} = \sigma_{yy} - \vec{\Sigma_{yx} \Sigma_{xx}^{-1} \Sigma_{yx}}^\intercal$ and $\vec{H} = \vec{X} - (\vec{\iota} \vec{\Gamma} + \vec{Z} \vec{\Delta})$. The marginal distribution of the treatment matrix, $ \vec{X} \sim MN( \vec{\iota \Gamma + Z \Delta }, \vec{I_n}, \vec{\Sigma_{xx}}) $ completes the joint distribution of $\vec{y}$ and $\vec{X}$.
This paper aims to incorporate model uncertainty into the framework outlined above (or its extension to non-Gaussian distributions as described in Subsection (ref)). In particular, a model refers to the exclusion of a specific set of covariates or instruments, or, equivalently, exact zero restrictions on the corresponding regression coefficients. We consider such uncertainty in both the treatment model and the outcome model.
Let $L \in \mathcal{L}$ denote a possible outcome model and let $M \in \mathcal{M}$ denote a possible treatment model, where $\mathcal{L}$ and $\mathcal{M}$ are the sets of all models considered. Then, the likelihood conditional on the models is given by
where the subscripts indicate the appropriate subset of the columns of the design matrices and equivalent zero restrictions on the coefficients. Note that the outcome also depends on the treatment model $M$ through the residual matrix $\vec{H}$. Both the number of potential outcome and treatment models are $2^p$. We do not introduce any cross-restrictions, so the number of combinations of outcome and treatment models is $2^{2p}$.
Unlike karl_instrumental_2012, we choose to always include the endogenous variables in the outcome model. This is a conscious choice and not necessary, as our prior specification and computational implementation would also easily allow for excluding the endogenous variables. We believe, however, that this is more natural as our primary focus is on estimating the effects of the endogenous variables while accounting for other covariates. This implies that we obtain a non-zero effect estimate in every model we consider. Compared to lee_incorporating_2022, our treatment model $M$ puts the same zero restrictions on all columns of the coefficient matrix $\vec{\Delta}$. This results in slightly less flexibility (for $l>1$) but leads to more interpretable results and substantially reduces the computational cost.
We also allow the user to fix certain variables to be instruments, meaning they are always excluded from the outcome model. In this variant, the user specifies two matrices: the pre-specified instruments $\vec{Z}$ and the covariates $\vec{W}$. The columns of $\vec{Z}$ are restricted to the treatment model, while the columns of $\vec{W}$ may enter either model. All of the results below continue to hold for this variant with slight modifications. Further details are provided in supplementary Section (ref).
The regression coefficient of $\vec{X}$ in the conditional model $\vec{y} \mid \vec{X}$ is generally not equal to $\vec{\tau}$. To see this, write the mean of the conditional distribution as
Naively regressing $\vec{y}$ on $\vec{X}$ targets $\left(\vec{\tau + \Sigma_{xx}^{-1} \Sigma_{yx}}^\intercal \right)$ instead of $\vec{\tau}$ itself. Whenever $\vec{\Sigma_{yx}} \neq 0$, this approach leads to biased results. This illustrates why it is necessary to consider the joint distribution of $\vec{y}$ and $\vec{X}$ instead of only the conditional distribution of $\vec{y} \mid {\vec{X}}$. Inference on $\vec\tau$ is still possible if we have suitable instruments (i.e. columns of $\vec{Z}_M$ that are not in $\vec{Z}_L$). Then, one can use the inference on $\vec{\Delta}$ from the treatment model to infer the covariance “ratio” $\vec{\Sigma_{xx}^{-1} \Sigma_{yx}}^\intercal$ and subsequently identify $\vec{\tau}$.
Our approach does not do this explicitly, but instead, in the outcome model, we condition on $\vec{H \Sigma_{xx}^{-1} \Sigma_{yx}}^\intercal$, which is known conditionally on the treatment and covariance parameters. The treatment residual $\vec{H}$ contains all the variation in $\vec{X}$ that the instruments and covariates do not explain. Therefore, including $\vec{H}$ acts as a control for the unobserved confounding. This is similar to classical control function approaches wooldridge_control_2015, but the coefficient of the control function is fixed (conditional on $\vec{\Sigma}$) and implied by our modeling assumptions. The treatment effect $\vec{\tau}$ is identified in a “classical” sense if at least $l$ instruments satisfy the following:
The exogeneity assumption above combines the unconfoundedness and exclusion restriction assumptions usually featured in a potential outcomes setup. The monotonicity assumption is trivially satisfied in our setup as the treatment equation is linear in the instruments imbens_instrumental_2014.
For a given combination of $L$ and $M$, the model is just-identified if the number of valid and relevant instruments implied by $L$ and $M$ is exactly equal to $l$. The model is under-identified if there are fewer instruments than $l$ and over-identified if there are more.
It is important to emphasize that identification is not a binary characteristic as it tends to be in classical inference for two reasons. First, if the model is under-identified, then the likelihood does not add information on some components of $\vec{\Sigma_{xx}^{-1} \Sigma_{yx}}^\intercal$. However, with a proper prior, we can still obtain a proper posterior distribution, albeit one that may be quite diffuse in some components. This carries through to the inference on $\vec{\tau}$. The second point is specific to our BMA approach: Each combination of outcome and treatment models can have different instruments with varying degrees of instrument strength. While some of these models will be well-identified, others might be very uninformative. Consequently, the marginal posterior of $\vec{\tau}$ can be very diffuse if the majority of the posterior weight is concentrated on these uninformative models. This merely reflects that we are not learning much about $\vec{\tau}$ from the available data, but it does not prevent us from conducting inference.
We can relax the Gaussianity assumption using the Univariate Link Latent Gaussian Models (ULLGM) framework proposed in steel_model_2024. The idea is to assign a latent Gaussian representation to any column of $\vec{[y : X]}$ that is not Gaussian and then perform posterior inference conditional on the latent Gaussian. Unlike lee_incorporating_2022, we do not need to rely on approximations.
More precisely, assume $\vec{y}$ and $\vec{X}$ are non-Gaussian and let the $n \times 1$ vector $\vec{q}$ and the $n \times l$ matrix $\vec{Q}$ be their latent Gaussian representations, respectively. Then, we model $(y_i, \vec{X_i}), \; i= 1,\ldots, n$, independently as
where
$F^{(y)}$ and $F^{(x_j)}$ are the distributions of $\vec{y}$ and $\vec{x_j}$, and $h_y$ and $h_{x_j}$ are invertible univariate link functions that map the latent Gaussian to the appropriate parameter of $F^{(y)}$ and $F^{(x_j)}, j=1\dots,l$. If required, $r_y$ and $r_{x_j}$ group any additional parameters of $F^{(y)}$ and $F^{(x_j)}$, which are assumed to be the same for all observations. We also assume the endogenous variables are independent conditional on their latent Gaussian representation, that is, all the dependence is in the Gaussian part. Many distributions can be expressed as members of the ULLGM family. Two examples that we use in our simulations and applications are steel_model_2024:
The interpretation of the Gaussian parameters varies with the distributions. For instance, in a Poisson-Log-Normal distribution for the outcome, the outcome regression parameters have a log-linear interpretation. The causal interpretation of $\vec{\tau}$ is then that of an expected log-ratio of expected potential outcomes (assuming $l=1$ for simplicity),
where $q_i(x) = \log \mathbb{E}_{y \mid q}\left[ y_i(x) \right]=\log \lambda_i(x)$ is the latent Gaussian potential outcome and $y_i(x)$ is the potential outcome on the observed level for a given value of $x$. Note that the causal interpretation of $\vec{\tau}$ is not affected by the choice of $F^{(x_j)}_{h_{x_j}(\cdot)}, j=1,\dots,l$, as $\vec{X}_i$ (and not $\vec{Q}_i$) multiplies $\vec{\tau}$ in the mean for $q_i$. If we choose to parameterize a location parameter of $\vec{X}_i$ by the latent Gaussian process $\vec{Q}_i$, then the invertibility of $h_{x_j}, j=1,\dots,l$, gives us a stochastic version of the monotonicity condition.
Model uncertainty is now introduced in the equations for the latent variables, similar to the discussion in Subsection (ref).
As a prior on the regression coefficients, we adopt a version of the widely used $g$-prior zellner_assessing_1986. However, we incorporate the intercepts into the $g$-prior as we cannot guarantee posterior propriety with an improper prior. In the outcome model, the prior distribution under model $L$ on the coefficient vector $\vec{\rho} = (\alpha, \vec{\tau}^\intercal, \vec{\beta}^\intercal)^\intercal$ then has the following continuous prior on $\vec{\rho}_L = (\alpha, \vec{\tau}^\intercal, \vec{\beta}_L^\intercal)^\intercal$, while the other elements of $\vec{\rho}$ are exactly zero:
where $\vec{U}_L = [\vec{\iota : X : Z}_L]$. In the treatment model, we use a matrix version of the $g$-prior on the coefficient matrix $\vec{\Lambda} = [\vec{\Gamma}^\intercal : \vec{\Delta}^\intercal]^\intercal$ given by a matrix normal prior on $\vec{\Lambda}_M = [\vec{\Gamma}^\intercal : \vec{\Delta}_M^\intercal]^\intercal$
where $\vec{V}_M = [\vec{\iota : Z}_M]$ and the other elements of $\vec{\Lambda}$ are zero. We will assume that for all models that we consider, the matrices $\vec{U}_L$ and $\vec{V}_M$ are of full column rank. This will ensure that the inverses in ((ref)) and ((ref)) exist for all models. Since the $g$-priors include the intercept, we recommend centering all the Gaussian components of $[\vec{y} : \vec{X}]$ to ensure that the zero mean prior on the intercept is reasonable. As stated in Subsection (ref), we impose the same zero restrictions on all treatment model equations if $l>1$. This means that the matrix normal prior in ((ref)) can be adopted, but in case we would want to allow potentially different covariates in each treatment model equation, we would require more flexible priors such as {\it e.g.} the Recursive Extended Natural Conjugate prior of RichardSteel. In case we would want to impose asymmetric shrinkage across treatment equations, we could adopt the Asymmetric Conjugate prior introduced by Chan22 in the context of VARs. Of course, the size of the model space and the complexity of the analysis would increase substantially and, for our purposes in this paper, we are typically not really interested in the coefficients in the treatment model. Hahn_etal_18 propose a factor shrinkage prior for the instrument coefficients when $l=1$, based on an approximation to the horseshoe prior (see Subsection (ref)).
The priors in ((ref)) and ((ref)) are conditionally conjugate for our sampling model, so we obtain closed-form expressions for conditional posteriors and marginal likelihoods (see Subsection (ref)). In addition, we only have to elicit two scalar hyperparameters, $g_L$ and $g_M$. The choice of $g_L,g_M>0$ controls the prior variance and the complexity penalty of the Bayes factors.
We can either fix $g_L$ and $g_M$ or put hyperpriors on them. We will focus on two possible choices: an adapted version of the benchmark or BRIC prior fernandez_benchmark_2001 and the hyper-$g/n$ prior liang_mixtures_2008. In the adapted benchmark prior, we fix $ g_L = \max \left\{n, (p+l+1)^2 \right\}$ and $g_M = \max \left\{n, (p+1)^2 \right\}$. In the context of the standard normal linear regression model, the benchmark prior results in model selection consistency fernandez_benchmark_2001, which will be examined in our setting in Subsection (ref).
The hyper-$g/n$ prior is characterised by its pdf, $$ p(g) = \frac{a-2}{2n} \left( 1 + \frac{g}{n} \right)^{-a/2}, $$ where $a>2$. This prior (unlike the regular hyper-$g$ prior) leads to consistent model selection in normal linear regression but does not yield analytic expressions for the marginal likelihood. Typical choices are $a=3$ and $a=4$, which behave similarly.
On any additional parameters $r_y$ and $r_{x_j}$ we adopt proper priors to ensure posterior propriety.
As in karl_instrumental_2012, we first put an inverse Wishart prior (as defined in supplementary Section (ref)) on the structural covariance matrix, $ \vec{\Sigma} \sim IW(\nu, \vec{I_{l+1}})$. Centering the prior over the (diagonal) identity matrix reflects the fact that we want inference on the degree of endogeneity to be primarily driven by the data. We recommend standardising the outcome and endogenous variables (or at least selecting appropriate units of measurement) to ensure the prior scale is not in strong conflict with the likelihood. The degrees of freedom parameter $\nu$ controls the amount of information in the prior and must satisfy $\nu > l$ so that the prior is proper. The choice of $\nu$ can be quite influential as {\it{i.a.}} it controls how tight the prior on all off-diagonal elements of $\vec{\Sigma}$ is around zero. It can substantially impact the inference on the marginal variances and the degree of endogeneity. To avoid being too dogmatic, we assume a hyperprior on $\nu$. An Exponential prior that is shifted by at least $l$, such that $\nu > l$, works well in our experience. For more numerical stability, it may be desirable to shift it slightly further. We use an Exponential with scale $1$ shifted by $l + 1$ as our default choice throughout the remainder of the paper. This prior is quite uninformative on $\nu$ itself while implying a prior on the covariance ratio $\vec{\Sigma_{xx}^{-1} \Sigma_{yx}}^\intercal$ that is (almost) identical to fixing $\nu = 3$ karl_instrumental_2012. Figure (ref) in the supplementary Subsection (ref) illustrates the implied priors on the covariance ratio and the outcome variance for shifted Exponential priors with different scales.
As a more flexible alternative, we also consider the Cholesky-based prior proposed in lopes_bayesian_2014. Write $\vec{\Sigma}$ as
where $\vec{a_{yx}} = \vec{\Sigma_{yx}} \vec{ \Sigma_{xx}}^{-1}$ denotes the scaled residual covariance. We specify the prior on $\vec{\Sigma}$ in three parts by assigning priors to the conditional outcome variance $\sigma_{y \mid x}$, the scaled covariance $\vec{a_{yx}}$, and the treatment variance $\vec{\Sigma_{xx}}$. This allows us to separately quantify the prior uncertainty about these individual covariance components, a potential advantage over the inverse Wishart prior used in karl_instrumental_2012.
A convenient choice is a Gaussian prior on the scaled covariance, $\vec{a_{yx}}^\intercal \sim N(0, \vec{\Omega}_a)$, an inverse Gamma prior (as in Section (ref)) on the outcome variance, $\sigma_{y \mid x} \sim IG(c, d)$, and an inverse Wishart prior on the treatment covariance matrix, $\vec{\Sigma_{xx}} \sim IW(\xi, \vec{I_l})$. These priors lead to conditional closed-form updates. As default choices for the variance hyperparameters, we use $c = \nu/2, d = 1/2, \xi = \nu-1$, where $\nu$ denotes the degrees-of-freedom parameter of the full inverse Wishart prior. Again, it can be useful to put a hyperprior on $\nu$. For the (scaled) residual covariance, we set $\vec{\Omega_a} = \omega_a \vec{I_l}$, where higher $\omega_a$ indicates a weaker prior belief in exogeneity. The key difference with the inverse Wishart prior is that we can now freely choose $\omega_a$, whereas it is restricted to $\omega_a = \sigma_{y \mid x}$ in the prior implied by the full inverse Wishart on $\vec\Sigma$.
We use a model prior based on the independent inclusion of variables in both the outcome and treatment equations. That is, we have inclusion probabilities, say $w_L$ and $w_M$, and the prior probability of model $L_j$ is $ p(L_j) = w_L^{p_j} (1 - w_L)^{p - p_j}$, where $p$ is the number of potential instruments and covariates and $p_j$ is the number included in model $L_j$. The choice of $w_L$ can be very influential, so we recommend a hyperprior on $w_L$ to make the procedure more adaptive. Using a Beta$(a,b_L)$ hyperprior results in the Beta-binomial prior
Following ley_effect_2009, we choose $a = 1$ and $b_L = (p-m_L)/m_L$, where the user specifies the prior mean outcome model size $m_L$. If not stated otherwise, we use $m_L = p/2$ as our default choice. For the treatment model we use a similar procedure, leading to uniform hyperpriors on the prior inclusion probabilities, as in ScottBerger10.
Note that this prior puts positive probability on all potential models. In particular, this means that we do not exclude underidentified models karl_instrumental_2012. As discussed above, in an underidentified model, the likelihood provides no information about certain components of the covariance. With a proper prior, however, we can still obtain a well-defined posterior, though some components may remain quite diffuse.
It also means that we allow for variables to enter the outcome model without entering the treatment model, which is generally discouraged in classical IV analysis angrist2009mostly. The reason is that the first-stage residuals are not necessarily uncorrelated with the covariates excluded from the treatment equation. This correlation leads to inconsistent coefficient estimates that can spill over into the treatment effect of interest. In our setting, however, this issue is less problematic for the following reason: If variables are excluded from the treatment equation, their correlation with the treatment is typically weak, meaning that any resulting second-stage bias has only a small effect on the treatment estimate. In addition, such a variable can still contribute to reducing residual variance in the outcome model and therefore increase precision.
A related question is what prior is induced on the number of valid and relevant instruments. This is described in supplementary Subsection (ref) and compared with the prior implied by IVBMA (with $l=1$).
Obtaining a tractable joint posterior distribution of all parameters and models across both the outcome and treatment models is not possible. Instead, we tackle the problem conditionally and take inspiration from the computational strategy in karl_instrumental_2012. We use a Gibbs sampler to iteratively sample the outcome parameters, treatment parameters, and the covariance matrix. Within the first and second steps, we update the models via a Metropolis-Hastings (MH) step in model space, update $g$ if it is not fixed, and then, given the model, we draw the parameters from their conditional posteriors. We describe each of these steps for updating the outcome model in more detail below, and the treatment model update is analogous. Algorithm (ref) summarises the Gibbs sampler.
\begingroup
\endgroup
Moving in model space is implemented through an MH step where the probability of accepting the proposed model $L'$ given the current model $L$ is
where $h$ is a proposal kernel. Throughout, models are proposed by randomly permuting the inclusion index of one covariate. This proposal is symmetric in the sense that $h(L \mid L')/h(L' \mid L) = 1$, so the acceptance probability reduces to the (conditional) Bayes factor times the prior ratio (capped at 1). This proposal can lead to slow mixing in very high dimensions, where other proposals might be more attractive, but we leave this to future work.
If $g_L$ is random, we update it using an MH step with a lognormal proposal. The proposal scale is tuned adaptively, targeting an acceptance rate of $0.234$.
Finally, we update the parameter vector $\vec{\rho} = (\alpha, \vec{\tau}^\intercal, \vec{\beta}^\intercal)^\intercal$ conditional on $L$ and $g_L$. The conditionally conjugate priors lead to a known distribution for the conditional posterior, as described in ((ref)) below in Subsection (ref).
We can obtain closed-form conditional posteriors and marginal likelihoods based on the prior specification described above. A detailed derivation is provided in supplementary Section (ref). In the outcome model, we have the conditional posterior
where $\vec{\Tilde{y}} = \vec{ y - H \Sigma_{xx}^{-1} \Sigma_{yx}}^\intercal$ is the endogeneity corrected outcome. The conditional Bayes factor (CBF) of model $L_i$ versus model $L_j$ is given by
where $\vec{U_i}, d_{U_i}$ and $\vec{U_j}, d_{U_j}$ are the design matrices and their number of columns in models $L_i$ and $L_j$ respectively.
In the treatment model, we obtain the (conditional) posterior
where $\vec{B_\Sigma} = \vec{I_l} + \frac{1}{\sigma_{y|x}} \vec{\Sigma_{yx}}^\intercal \vec{\Sigma_{yx} \Sigma_{xx}^{-1}}$ and $\vec{\Tilde{X}} = \vec{X} - \frac{1}{\sigma_{y|x}} \vec{\epsilon \Sigma_{yx} }\left(\vec{B_\Sigma^{-1}}\right)^\intercal$. The CBF of model $M_i$ versus model $M_j$ is
where $\vec{A_\Sigma} = \left( \left( \vec{I_l} + g_M^{-1} \vec{B_\Sigma}^{-1} \right)^{-1} \right)^\intercal \vec{\Sigma_{xx}}^{-1} \vec{B_\Sigma}$ and $\vec{V}_i, d_{V_i}$ and $\vec{V}_j, d_{V_j}$ are the design matrices and their number of columns in models $M_i$ and $M_j$ respectively.
The treatment model uses the joint distribution of $\vec{y}$ and $\vec{X}$ for its marginal likelihood, while the outcome model uses the conditional distribution $\vec{y} \mid \vec{X}$. This distinction arises because the marginal distribution of $\vec{X}$ is independent of outcome parameters, but the outcome model $\vec{y} \mid \vec{X}$ depends on the treatment parameters.
Using the IW prior the conditional posterior for the covariance matrix is
Conditional on the data and the outcome and treatment parameters, the residuals $\vec{\epsilon}$ and $\vec{H}$ are known, so the posterior scale matrix can be computed within the Gibbs steps.
For the Cholesky-based prior, we can use the conditional posteriors for the three covariance components,
A draw of the covariance matrix $\Sigma$ is then implied by equation ((ref)).
If we put a hyperprior on $\nu$, we add an extra step to the Gibbs sampler updating $\nu$. Given $\vec{\Sigma}$, $\nu$ is independent of everything else, so the full conditional is proportional to the prior on $\vec{\Sigma}$ given $\nu$ times the prior on $\nu$. We use an adaptive MH step, targeting an acceptance rate of $0.234$.
Here, we discuss the computational strategy to deal with the ULLGM models introduced in Subsection (ref). We add an extra step to the Gibbs sampler that draws the latent Gaussian (and potentially additional parameters $r_y$ or $r_x$) and then perform posterior inference on the Gaussian parameters conditional on the latent Gaussian representation. steel_model_2024 propose to use an MH step with an adaptive Barker proposal livingstone_barker_2022 to sample the latent Gaussians. This is a good compromise between using gradient information to increase the mixing speed and maintaining robustness.
Define the residuals based on the latent Gaussian $\epsilon_i = q_i - \left(\alpha + \vec{X_i \tau} + \vec{Z_i \beta} \right)$ and $\vec{H_i} = \vec{Q_i} - \vec{\Gamma}^\intercal - \vec{Z_i \Delta}$. For the outcome, the gradient used in the Barker proposal is
where the first term depends on the specific distribution of $y_i$, and the second term comes from the Gaussian “prior”. For the $j$-th endogenous variable, the gradient is
The first term depends on the distribution of the $j$-th endogenous variable, the second arises from the outcome distribution, and the last is the contribution of the Gaussian “prior”.
If additional parameters $r_y$ or $r_{x_j}$ are required, we add an extra step to update them separately after the respective component of the latent Gaussian. This is done through an adaptive MH step, targeting an acceptance rate of $0.234$.
After updating the latent Gaussian representation and any additional parameters, we can proceed with the same steps as in the Gaussian case. The formulas from the latter remain valid using the latent Gaussian representation in place of any non-Gaussian variables.
It is important to know that if the sample size $n$ increases, we end up putting more and more posterior mass on the “correct” model. We are assuming here that the model that actually generated the data lies within the model space considered.\footnote{This is what is often referred to as an $M$-closed setting. In most situations, model selection consistency naturally extends to the $M$-open framework by selecting the model that is closest to the true model according to a suitable metric Mukhopadhyay. We leave a detailed treatment of this for future work. }
Let us assume that the data are generated from models $L_i$ and $M_i$ for, respectively, the outcome and treatment equations. Then we say that our gIVBMA procedure is model selection consistent if both CBF($L_i,L_j$) and CBF($M_i,M_j$) tend to $\infty$ with $n$ for any $L_j\ne L_i$ and $M_j\ne M_i$. We show that the following results hold:
For both the fixed (BRIC) and random $g$ (hyper-$g/n$) options used throughout this paper, we can thus show that model selection consistency holds. This result provides additional justification for data-driven instrument selection. If the sample size is large enough, posterior mass concentrates on the “correct” models and, therefore, instruments and covariates will be correctly separated.
We evaluate the performance of gIVBMA in three different scenarios and compare it to naive BMA, IVBMA, and several classical methods. Throughout, we consider the median absolute error (MAE) and median bias of the point estimates of $\vec{\tau}$, (credible or confidence) interval coverage of $\vec{\tau}$, and the log predictive score (LPS) on a separate holdout dataset. For all Bayesian methods, we use the posterior mean as our point estimate for the MAE and bias calculations. While MAE, bias, and coverage relate specifically to the quality of the estimation of $\vec{\tau}$, LPS is a measure of predictive adequacy, focused on probabilistic predictions of the outcomes. In each scenario, we simulate $100$ independent datasets of size $n$, and for the LPS computation, we also generate holdout datasets of size $n/5$. Supplementary Section (ref) provides more details on the different measures used, while Section (ref) describes the different estimators and their implementation.
The predictive comparison is perhaps not entirely fair to some of the classical methods (in particular those that are two-stage in nature) as they reduce the bias in $\vec{\tau}$ but lack a correction term like gIVBMA. Consequently, when making predictions, these methods use an estimate of $\vec{\tau}$ instead of $\vec{\tau + \Sigma_{xx}^{-1} \Sigma_{yx}}^\intercal$, which would be the appropriate coefficient for the conditional model $\vec{y} \mid \vec{X}$. This can lead to poor predictions when $\vec{\Sigma_{xx}^{-1} \Sigma_{yx}}^\intercal$ becomes large. However, in practice, these methods might not be employed when prediction is the primary objective. Methods that do not correct for endogeneity (e.g., naive BMA) target the biased coefficient $\vec{\tau + \Sigma_{xx}^{-1} \Sigma_{yx}}^\intercal$, but provide good predictions as this is the right coefficient for predicting $\vec{y}$ conditional on $\vec{X}$. This discrepancy reflects the fact that predicting observed and counterfactual outcomes are different problems under endogeneity. Our probabilistic method allows us to provide endogeneity-corrected inference on $\vec{\tau}$ while maintaining good predictive power in the observational model.
To investigate the performance in the presence of invalid instruments, we consider a simulation setup similar to the ones in kang_instrumental_2016 and windmeijer_use_2019. We consider two different sample sizes $n \in \{50, 500\}$ and $p = 10$ potential instruments simulated from independent standard Normals. Their coefficients in the outcome model are set to $\vec{\beta} = (1, \ldots, 1, 0, \ldots, 0)$ where the first $s$ coefficients are one (and, therefore, the first $s$ instruments are invalid). We vary $s \in \{3, 6\}$ such that the plurality rule only holds in the first scenario. The instruments' treatment coefficient is set to $\vec{\delta} = c \vec{\iota}$, where $c>0$ is chosen such that the first-stage $R^2_f$ is approximately $0.2$. The treatment effect is set to $\tau = 0.1$. Then the data is generated from ((ref)), while $\sigma_{yy}=\sigma_{xx}=1$ and $\sigma_{yx}=1/2$.
To compare the performance of our different prior specifications, we evaluate several variants of our method using BRIC or hyper-$g/n$ priors on the regression coefficients, and either the full inverse-Wishart or the Cholesky-based covariance prior, the latter with different hyperparameters $\omega_a \in \{0.1, 1, 10\}$. We benchmark our approach against IVBMA karl_instrumental_2012, TSLS with all instruments, oracle TSLS, and the sisVIVE estimator kang_instrumental_2016. We also use a Bayesian IV method that includes all instruments in both models with horseshoe priors carvalho2009handling on their regression coefficients. See supplementary Section (ref) for more details.
Table (ref) reports the instrument selection performance of the different gIVBMA variants. The true number of valid and relevant instruments implied by the data-generating process is \(N_Z = p - s\). For \(n = 50\), all variants distribute their posterior mass between \(N_Z = 0\) (exceeding the prior probability) and \(0 < N_Z < p - s\), placing only negligible mass on the true value (or above). When \(s = 3\), the mass is roughly evenly split, whereas for \(s = 6\), the models assign more weight to non-identified specifications. In contrast, for \(n = 500\), all variants allocate only negligible posterior mass to non-identified models. Overall, the hyper-\(g/n\) specifications yield consistently better performance, while smaller values of $\omega_a$ seem to do better for $n=500$.
Table (ref) shows that all gIVBMA variants predict well and lead to small bias and MAE. For small $n$, the BRIC specifications mostly have a higher MAE and median bias than the hyper-$g/n$ variants. The IVBMA method also predicts well and has low MAE and median bias, but overcovers for small $n$ and significantly undercovers for large $n$. TSLS results in a very large positive bias, very small coverage, and poor prediction. The oracle version of TSLS does better but is slightly worse than the gIVBMA methods for small $n$. All the methods above appear largely unaffected by the plurality rule. The latter does affect sisVIVE, however, which performs very poorly when $s=6$, as expected. The Bayesian horseshoe (BayesHS) method incurs substantial bias and does not cover well at all.
In the subsequent simulations, we shall focus on the hyper-$g/n$ prior on the regression coefficients, using either the full inverse-Wishart covariance prior or the Cholesky-based prior with $\omega_a = 0.1$, as these performed best in this experiment. The Cholesky prior with $ \omega_a = 1$ also performed well, but behaved similarly to the full inverse-Wishart prior. This is not surprising, since $\sigma_{y \mid x}$ is roughly of the order 1 in the IW prior (see Supplementary Figure (ref)).
We consider an example with two endogenous variables, a Gaussian and a Beta, and correlated instruments. We generate 15 valid instruments (five of which are also relevant) with a correlation structure similar to fernandez_benchmark_2001. The data-generating process is given in detail in supplementary Section (ref). We vary the sample size $n \in \{50, 500\}$.
We compare gIVBMA against naive BMA, IVBMA, OLS, TSLS (using all instruments and an oracle version), and MATSLS. The latter refers to the model-averaged TSLS estimator with unrestricted weights by kuersteiner_constructing_2010 and is included as it is the only feasible classical method that is competitive with oracle TSLS in all scenarios of Subsection (ref). IVBMA cannot directly model the Beta-distributed endogenous variable and instead relies on a Gaussian approximation.
Supplementary Table (ref) presents the full estimation error, coverage, and LPS results. Both retained gIVBMA variants perform best in terms of MAE and median bias, for both sample sizes. BMA is severely biased as expected, but coverage for the second endogenous variable is adequate (in contrast with the coverage for $\vec{X}_1$). This is a consequence of the non-linear transformation from the latent Gaussian to the Beta variable, diluting the endogeneity. TSLS, oracle TSLS, and MATSLS all perform similarly well. IVBMA results in a substantial bias for both sample sizes, yet it predicts and covers relatively well.
Figure (ref) presents the distributions of the posterior probability of the true treatment model and the mean treatment model size. Unlike gIVBMA, IVBMA allows for separate treatment models for the two endogenous variables (see Subsection (ref)), and we report the results separately for $X_1$ and $X_2$. Both gIVBMA variants put substantial posterior mass on the true model, and their mean model size is accordingly very close to that of the true model (which is $5$). For $n=50$, IVBMA tends to select too many variables for $X_1$ and too few for $X_2$. The performance of IVBMA is better in the $n=500$ setting, but it still cannot compete with gIVBMA. Supplementary Table (ref) presents median posterior inclusion probabilities (PIP) for all instruments. The selection performance of gIVBMA is clearly superior for both sample sizes.
We also consider a setting with many weak instruments, similar to the simulation study in kuersteiner_constructing_2010, and a variant with a Poisson endogenous variable. In this case, the instruments are assumed to be valid a priori and cannot be included in the outcome model, using the gIVBMA version described in Section (ref). In the Gaussian setting, the gIVBMA variants deliver good coverage and predictions, but are worse than IVBMA in terms of MAE. By contrast, in the Poisson setting, gIVBMA is clearly superior in MAE, and IVBMA has very poor coverage. The horseshoe prior yields lower MAE than gIVBMA in most scenarios, but tends to overcover and predicts poorly. In addition, gIVBMA (like all the other methods with exact zero restrictions) formally deals with identification and has an easier interpretation. MATSLS is competitive for large $n$ but predicts substantially worse than gIVBMA throughout. Full details are provided in supplementary Section (ref).
carstensen_primacy_2006 show that geographic factors such as disease ecology, particularly malaria prevalence, have a substantial negative impact on income. ditraglia_using_2016 reanalyzes these data to illustrate the usefulness of their proposed Focused Moment Selection Criterion. They consider the regression model $ \log \text{gdpc}_i = \beta_1 + \beta_2 \text{rule}_i + \beta_3 \text{malfal}_i + \epsilon_i$, where $\text{gdpc}$ is real GDP per capita in 1995 prices, $\text{rule}$ is an average governance indicator measuring the quality of institutions, and $\text{malfal}$ is the fraction of the population at risk of malaria transmission in 1994. The quality of institutions and the prevalence of malaria are likely endogenous. The dataset (for 44 countries) contains various potential instruments such as historical settler mortality (lnmort), malaria transmission stability (maleco), winter frost levels (frost), the maximum temperature during peak humidity (humid), latitude (distance from equator), the proportion of Western Europeans and English speakers (eurfrac and engfrac), proximity to coast (coast), and predicted trade share (trade). Settler mortality (lnmort) and malaria ecology (maleco) are very plausible to be exogenous, but there is uncertainty about the exogeneity of all other potential instruments.
We will consider both institutional quality and malaria prevalence as joint endogenous variables (i.e. $l = 2$) and compare the PIP of all instruments in the outcome and treatment model to assess their validity. The rule variable is approximately Gaussian. The malfal variable is a proportion and only takes values in $[0, 1]$, so we model it using a Beta distribution. We put an Exponential prior with rate $1$ on the additional dispersion parameter of the Beta distribution and draw it in an MH step with a Gaussian proposal. One challenge is that several countries in our sample have recorded malfal values of exactly 0 or 1. This is a problem for our algorithm as the logistic function (e.g., used in computing the gradient for the proposal) is not finite at 0 or 1. We deal with this incompatibility of our observations with the sampling model by using set observations fernandez_multivariate_1999. We treat zero observations (and analogously ones) as belonging to a set $(0, 0.0005)$ and add an extra step to draw the value from its sampling distribution truncated to that set. This matches these observations to the dominating measure (Lebesgue) of the sampling model.
Table (ref) shows the results. All outcome PIPs are very low, while four of the instruments are included in the treatment model with high probability, indicating strong support for overidentification. In particular, maleco and lnmort are valid and relevant instruments, as expected, while the humidity and coast variables are valid and highly relevant as well. The hyper-$g/n$ prior tends to yield a larger $g$ in the outcome model and a smaller $g$ in the treatment model than the fixed $g$ corresponding to BRIC, leading to lower and higher PIP, respectively. Our treatment effect estimates are in line with those in ditraglia_using_2016. Their point estimates for rule and malfal are between $0.81$ and $0.97$ and between $-1.16$ and $-0.9$, respectively, depending on which instrument set they use. The posterior distributions of both components of the covariance vector, $\sigma_{12}$ and $\sigma_{13}$, concentrate near zero. Thus, the gIVBMA results are similar to naive BMA, though gIVBMA credible intervals are slightly narrower, reflecting the borrowing of strength from the treatment equation.
To assess the predictive performance on this dataset, we compute the LPS using leave-one-out cross-validation. For each observation, the model is trained excluding that observation, and the LPS is then calculated for the excluded data point. We compare the performance of gIVBMA, IVBMA, BMA, and TSLS. Although sisVIVE would be relevant given potential concerns about instrument validity (at least a priori), we exclude it as this approach does not accommodate multiple endogenous variables. Methods designed for many weak instruments were also omitted since their predictive performance was consistently inferior to TSLS in our simulations. Table (ref) presents the mean LPS values. The gIVBMA variants perform best, closely followed by naive BMA and TSLS, while IVBMA predicts considerably worse.
A prominent problem in microeconomics is estimating the causal relationship between educational attainment and earnings. Education levels are not randomly assigned but are likely influenced by unobservable characteristics that simultaneously affect earning potential. Therefore, the observed correlation between education and wages may not accurately reflect the true economic returns to schooling. Many econometricians have tried to isolate the causal effect by finding instrumental variables that provide plausibly exogenous variation in educational attainment, such as quarter of birth angrist_does_1991 or geographic variation in college proximity card1995collegeproximity.
We revisit this problem using the dataset from card1995collegeproximity, a subset of the National Longitudinal Survey of Young Men, containing all men who were interviewed in 1976 and provided valid wage and education responses. Our outcome is the logarithm of hourly wages with years of schooling as the single endogenous explanatory variable of interest. Our set of potential exogenous covariates and instruments includes experience, experience squared, college proximity (distance to a 2 or 4-year college), variables on family background, marital status, race, and regional indicators. The dataset also includes information on parents' educational attainment, which have a substantial proportion of missing values. As in card1995collegeproximity, we impute those missing values using their mean and include indicators for missingness. After imputing missing values, the data consists of $3,003$ observations and $23$ potential instruments and exogenous controls (see supplementary Table (ref)). We use fully data-driven instrument selection such that all variables can be used as instruments or covariates.
Figure (ref) presents our posterior results.\footnote{Fitting each specification requires less than an hour on a standard 12th Gen Intel i7 CPU. Thus, computational costs are not a limiting factor for this approach, making gIVBMA practical for datasets with several thousand observations.} Consistent with card1995collegeproximity, the gIVBMA posterior for the treatment effect $\tau$ indicates higher returns to schooling than BMA models that lack an endogeneity correction. For comparison, the TSLS estimate when using only proximity to a four-year college as an instrument and all other variables (except proximity to a two-year college) as exogenous controls is $\hat{\tau} = 0.142$ ($95\%$ CI: $[0.028, 0.256]$). Thus, our posterior results support slightly lower returns to schooling with higher estimation precision, a discrepancy likely driven in part by differences in instrument selection. While our algorithm frequently selects proximity to a four-year college as an instrument, it always includes parental education as well. In contrast, IVBMA yields implausible results by concentrating most of its posterior mass exactly at zero. This might occur because IVBMA fails to select experience in the outcome model and uses it as an instrument instead (see supplementary Table (ref)). IVBMA tends to have a higher number of instruments as it selects fewer variables in the outcome model (see supplementary Figure (ref)). For all methods, the posterior probability of non-identification (no instruments) is zero.
As a robustness check, we re-estimated the model using only observations with complete parental education data ($n=2,215$; see supplementary Section (ref)). The results remain qualitatively similar to the main specification. The primary difference is that college proximity is selected as an instrument less frequently, while parental education remains a constant choice. Notably, the IVBMA results appear more plausible in this subsample, as the model shifts toward using parental education as its primary instrument (see supplementary Table (ref)).
Supplementary Table (ref) shows that the predictive performance of gIVBMA, IVBMA and naive BMA is relatively similar on both versions of the dataset, while TSLS does worse.
We have introduced the gIVBMA method, which performs posterior inference in structural equation models with at least one endogenous regressor, averaging over all possible sets of potential instruments and covariates. We allow for the model to freely choose its instruments based on the data, which provides additional robustness against invalid instruments. Our computational strategy relies on updating the outcome model, the treatment model, and the covariance matrix separately using conditional Bayes factors in the model updates. Exploiting the ULLGM framework steel_model_2024, we can accommodate non-Gaussian outcomes and endogenous variables. We provide simple and unrestrictive necessary and sufficient conditions for model selection consistency of gIVBMA. In both simulation and real data experiments, gIVBMA outperforms the earlier method proposed by karl_instrumental_2012 and many classical estimators considered in the literature.
Future work could extend the framework to more complex covariance structures, possibly rendering gIVBMA available for settings such as time series or network data.
All derivations and proofs, along with further details on the variant with fixed instruments, prior specification, simulation experiments, and empirical examples, are provided in the Supplementary Material. The code and data used are available at https://github.com/gregorsteiner/gIVBMA-Code.