Applying Bayesian Multivariable Mendelian Randomisation to Prioritise Candidate Causal Traits From High‐Dimensional Data: Illustration From Estimation of the Effect of Maternal Metabolites on Offspring Birthweight

1.1. MVMR‐BMA

In this section, we explain MVMR‐BMA, the inputs required, and the interpretation of the outputs with reference to our application throughout. Figure 1 gives the summary of suggested steps when running MVMR‐BMA.

Bayesian model averaging is a methodological approach to estimation which incorporates uncertainty in the model specification, thus reflecting that the true causal model is a priori unknown. Within a Bayesian model averaging approach an algorithm searches over a range of different models to estimate the probability with which each is the likely true model and incorporates model uncertainty in the effect estimate. The different models considered could be different exposures, but could also be various specifications of other aspects of the model, such as whether a relationship is linear or whether a non‐linear model should be used. It is a hypothesis‐free approach that avoids the need for the researcher to specify a single assumed ‘true’ model prior to estimation.

MVMR‐BMA uses a Bayesian model averaging approach to select which exposures out of a larger set of, potentially correlated, exposures are causal exposures for the outcome, within an MR framework (Zuber, Colijn, et al. ²⁰²⁰). For each model, the effects are estimated using conventional Bayesian MVMR approaches based on conjugate Normal priors and thus are in spirit close to Ridge regression. The models are estimated using a pre‐specified set of SNPs each of which are associated with at least one of the potential exposures, using the same specifications for association as would be applied in a `conventional MR estimation (Sanderson et al. 2022, ²⁰¹⁹; Burgess and Thompson ²⁰¹⁵). MVMR‐BMA considers all possible combinations of exposures from the initial set (up to a pre‐defined maximum included in any one model) and incorporates model uncertainty through calculating a weighted average of the effect estimate from each model (Hoeting et al. 1999). Every model considered is given a model‐specific posterior probability (PP) of being the true causal model. The PP is the probability of how likely the model is to be the true causal model, enabling comparison of the different models considered by the method. The PP is estimated based on a combination of prior probability for that particular model and the Bayes factor calculated from the data. Bayes factor is a commonly used statistic for Bayesian model comparison, which compares the likelihood of one model relative to another (Zuber, Colijn, et al. ²⁰²⁰).

The BMA algorithm runs for many iterations, selecting a different model each time but returning to some models more often due to a relatively high likelihood of being the true model calculated by the estimation (and/or high prior probability specified by the researcher of being the true model). While for a smaller number of exposures, the traditional BMA (Hoeting et al. 1999) can exhaustively evaluate all models, that is, all possible combinations of exposures, this becomes computationally infeasible with standard desktop computing resources when the number of exposures approximately exceeds 20 (total number of possible models is 2²⁰). Zuber et al. (Zuber, Colijn, et al. ²⁰²⁰). have used shotgun stochastic search to reduce the burden. This allows the method to be applied to sets of exposures that are much larger than 20 without excessive computational burden. The Supplementary Notes and Supporting Information S1: Figure S1 compares runtime for different algorithm using simulated data.

For each individual exposure, MVMR‐BMA aggregates the evidence across models in which the exposure is included and calculates the marginal‐inclusion probability (MIP) and the model‐average direct causal effect. The MIP for each exposure is an estimate of how likely that exposure is to have a causal effect on the outcome and is calculated as the sum of PP across all models that contain the corresponding exposure. The model‐averaged direct causal effect of the exposure on the outcome is the effect estimate for each exposure, averaged across the models that include the corresponding exposure (Zuber, Colijn, et al. ²⁰²⁰).

1.2. Model Input Considerations

1.3. Instrument Selection

In MR, genetic variants for exposures are selected on the basis that they have a statistically strong association with the exposure(s) of interest, usually defined as being genome‐wide statistically significantly associated (p < 5 × 10⁻⁸) with the exposure in a GWAS of the exposure, pruned so only uncorrelated (R² for linkage disequilibrium < 0.001) variants are retained (Sanderson et al. 2022). In MVMR, this is extended to select all genetic variants that are strongly associated with any of the exposures using the same criteria for selection as for MR (Sanderson et al. 2019).

We applied this approach to the selection of genetic variants as instruments in our estimation. We selected all variants associated with our final set of 46 traits with a p‐value < 5 × 10⁻⁸ in a GWAS of the full set of metabolites that had previously been conducted in UK Biobank (Borges et al. 2022) (summary statistics available in the IEU OpenGWAS (Bycroft et al. 2018; Elsworth et al. ²⁰²⁰; Hemani et al. ²⁰¹⁸). This list was pruned to remove variants with an R ² for linkage disequilibrium < 0.001 across the full list of variants associated with any exposure.

1.4. Outcome Data

For our outcome data set, we used a GWAS meta‐analysis including 210,267 European participants, combining results of 41 studies from the EGG (Early Growth Genetics) consortium and UKBB (UK Biobank), with a mean age at delivery ranging from 24.5 to 31.5 years (Warrington et al. 2019). This GWAS was adjusted for foetal genotype to avoid potential biases due to correlated with the maternal genotype (Beaumont et al. 2018). From this data set, we extracted the association between each of the SNPs identified above and offspring genotype to include in our estimation.

1.5. Priors

Like all Bayesian analysis, MVMR‐BMA requires prior parameters to be defined, which are probability and variance. Prior probability in this model defines the probability that any individual exposure will be a causal exposure. A smaller value informs the algorithm to favour models with fewer exposures, and a larger value to favour models with many exposures. Zuber et al. (Zuber, Colijn, et al. ²⁰²⁰) recommend specifying the probability based on the expected model size (i.e., probability times number of exposure). The causal effect estimate follows a normal prior centred at zero with a user defined variance of σ ². This specification informs the method whether to favour exposures with larger or smaller causal effect estimates, with a larger σ ² favouring exposures with larger (i.e. further away from the null) causal effects.

In the main analyses, we ran MVMR‐BMA using a prior probability of 0.1 for each of our 46 exposures being a causal exposure for birthweight, as we are expecting that the model will contain four to five metabolites (= 0.1 $\times$ 46). The prior variance was set to 0.5 for each metabolite, as we would like to allow for stronger causal effects. We also ran sensitivity analyses of different choices of prior probability and variance.

3.1. Priors

It is possible that the chosen prior will have a substantial influence on the result obtained from the MVMR‐BMA, and therefore the estimation should be run for a range of priors. If a model is consistently selected across priors, it can be considered to be robust as the prior itself is not determining the selected model. We therefore tested a range of values for prior probability and prior variance for our model. All the results are in Supporting Information S1: Table S4 and S5. For prior probability, we considered values of 0.01, 0.05 and 0.20, and for prior variance, we considered values of 0.30, 0.50 and 0.70. As expected, with a smaller prior probability, the model with fewer exposures (only glucose) is favoured over the model with more exposures (glutamine and glucose). When a larger prior probability is used, the model with both glutamine and glucose is favoured over just glucose (Supporting Information S1: Table S4). A prior variance of 0.7 led the BMA to select larger values for the causal effect estimate, in our case glucose has the largest effect (model‐specific causal effect of 0.24 for the model only including glucose) compared to other metabolites, and therefore the model with glucose only becomes most favoured by BMA PP changed from 0.311 to 0.41 (Supporting Information S1: Table S5a). However, the model with glutamine and glucose and the model with glucose are consistently favoured from the MVMR‐BMA procedure. Therefore our BMA analyses are robust to differing values of prior variance and prior probability.

3.2. Convergence and Mixing

Initially our MVMR‐BMA algorithm was run using 100,000 iterations throughout, however we additionally re‐ran the algorithm with 500,000 iterations to check that the algorithm was being run for sufficient iterations. We obtained identical results to the main analysis (Supporting Information S1: Table S6). We conclude 100,000 iterations were sufficient for convergence and mixing. Thus we used 100,000 iterations for all of the sensitivity analyses.

4.1. Outliers and Influential Variants

MR estimation relies on a set of IV assumptions for the estimates obtained to be interpreted as causal effects of the exposure(s) on the outcome (Sanderson et al. 2022). One assumption is that the SNPs used as instruments only affect the outcome through the exposure and not by any other pathway. If the majority of the SNPs used in the estimation are valid (in the sense that they satisfy this condition), then SNPs that violate the assumption will give individual effect estimates that differ from those for the larger set of SNPs. These SNPs will be outliers if the individual SNP level estimates are considered, and will make the total heterogeneity in the effect estimates across all SNPs larger. One way heterogeneity in the model can be measured is through a Q‐statistic. The appropriate critical value for this test can be obtained from a ${\chi }^2$ distribution with degrees of freedom equal to the number of SNPs minus the number of exposures. A statistically significant value for the Q‐statistic (i.e., larger than the relevant critical value 0) indicates heterogeneity in the model is higher than would be expected purely by statistical noise (Bowden et al. 2018).

Similarly, particularly influential SNPs that have more influence on the results obtained, may make the results from the analysis less robust if those SNPs are the ones that violate the IV assumptions. Cook's distance is a measure of the influence of any individual data point (in this case, SNP, as we are conducting analysis with summary statistics) on the estimates obtained from a regression analysis. In this setting, each SNP is a data point, and the IVW estimation is the regression in question.

To assess the influence of outliers and influential SNPs in MVMR‐BMA we re‐ran the analysis excluding any SNPs that were found to be outliers or highly influential in the top model (including glutamine and glucose). SNPs were defined as outliers where the individual SNP‐level contribution to the overall Q statistic was greater than the relevant critical value (here 22.78). Highly influential SNPs were defined to be those exceeding the Cook's distance threshold for the IVW estimates (here 0.987). Note: these thresholds are dependent on the number of SNPs and so will differ across applications. A total of three outlying and zero highly influential SNPs were removed based on these criteria. Q‐statistics and Cook's distance are reported in Supporting Information S1: Tables S7 and S8 and Figure S2.

This exclusion changed the top models selected, the model with the highest PP was now the glucose‐only model, while the glucose and lactate model ranked second, which was previously not in the top five. Correspondingly, the MIP for lactate increased, and the MIP for glutamine decreased. The SNPs removed were rs1801133, rs2168101 and rs7137828. Outlying SNPs should be removed from any MR analysis with caution. A SNP being outlying does not necessarily mean that it is pleiotropic, and it is possible that outlying SNPs are the SNPs that give the true causal effect. Investigation of the identified outlying SNPs showed that rs2168101 is strongly associated (p < 5 × 10⁻⁸) with 5 of the 46 metabolites, including both glutamine (p = 7.6 × 10⁻⁴⁸) and alanine (p = 9.0 × 10⁻³¹). rs1801133 is only associated with glycine (p = 1.8 × 10⁻⁸). rs7137828 is strongly associated with 15 traits (p < 5 × 10⁻⁸) including lactate (p = 6.8 × 10⁻²⁴). It is not possible to determine from this data alone whether the removal of these SNPs has removed pleiotropic effects in the estimation or if the removal of these SNPs has changed the results obtained because the ability of the SNPs used as instruments to predict glutamine, alanine and lactate has changed between the two models. However, the result that glucose is the top exposure in each estimation adds strength to the conclusion that glucose is an important determinant of birthweight based on genetic evidence.

4.2. Exposure Inclusion

Two exposures, total lipids in very large VLDL and total lipids in small LDL, were excluded from the main analysis due to their high genetic correlation with total lipids in large VLDL and total lipids in medium LDL, respectively. Therefore, as sensitivity, we excluded the total lipids in large VLDL and total lipids in medium LDL instead, to include their correlated counterpart. Similar results to the main analysis were obtained (Supporting Information S1: Table S11).

5.1. Instrument Strength in MVMR‐BMA

A key assumption for MR and MVMR estimation is that the exposure(s) can be strongly predicted by the genetic variants used as instruments in the relevant population for the estimation. F‐statistics in for the association between the instruments and the exposure can be used to test the strength of the instruments in an estimation (Sanderson et al. 2022). A value of 10 or higher is often used as an indicator that there is unlikely to be substantial weak instrument bias in the estimation, based on an approximation of the critical value for a standard test for weak instruments (Stock and Yogo 2005). This relies on the assumption that the data where the exposure is measured in is the relevant population for the outcome, and so in the example here, assumes that the association between the genetic variants and the metabolites of interest is the same in pregnant women (the relevant population for outcome data) as it is in the general population (in which our exposure data was measured).

MVMR‐BMA includes all SNPs associated with any candidate trait in each model considered. If only very few of the SNPs included are associated with a trait that has a causal effect on the outcome, then most of these SNPs will show evidence of no causal effect on the outcome. This may lead the MVMR‐BMA algorithm to struggle to distinguish between noise, driven by the many SNPs associated with non‐causal traits, and the true effect, which is seen only for a few SNPs. In our application, 17 of the 598 SNPs included in our main analysis are associated with glucose at genome‐wide significance level. For lactate, which was only selected by our sensitivity analyses, 7 of the 598 SNPs were strongly associated (Supporting Information S1: Table S9). Weak instruments in MVMR estimation can bias the causal effects either towards or away from the null, however, they also lead to higher uncertainty around the estimates obtained. This can lead to lower support for a model being the true causal model, down weighting the PP given to the model in a MVMR‐BMA approach.

The issue of weak instruments is not specific to MVMR‐BMA but will appear, in slightly different forms, in any MVMR analyses attempting to estimate multiple traits when only a few SNPs are associated with an exposure of interest (Sanderson et al. 2021). We calculated conditional F‐statistics, see Supporting materials ‘Conditional F‐statistics’, for each of the top models selected by the MVMR‐BMA to assess the influence of weak instruments on model selection (Sanderson et al. 2021). Conditional F‐statistics for each exposure ranged from 8.3 to 83.6 (Supporting Information S1: Table S10) when only SNPs associated with the selected exposures were considered. However, when all SNPs are included in the estimation, as is the case in MVMR‐BMA, these ranges dropped to 2.6 to 12.5. This difference occurred because the F‐statistic is calculated as an average across all genetic variants used as instrumental variables in the estimation. Glucose had a conditional F‐statistic of approximately 4.6–4.9, depending on exact model specification, and so the MVMR‐BMA effect estimates obtained are likely to have been biased by weak instruments. However, due to the set up of MVMR‐BMA, the causal effects are biased towards the null. Notably, lactate, which was not identified as having a high MIP in the primary analysis, had a conditional F‐statistic of approximately 2.7 in the models in which it was included. Effect estimates for both lactate and glucose are therefore likely biased by weak instruments, thus caution should be taken when interpreting their effect estimates (Sanderson et al. 2021). Checking the instrument strength of the exposures in the top models is important in interpreting the results from any MVMR‐BMA estimation. This issue applies similarly across MVMR estimation and currently there are limited approaches available for weak instrument robust MVMR estimation (Sanderson et al. 2021).

5.2. What to Look for to Check the Model Is Working

The MVMR‐BMA algorithm will always rank the models and exposures based on their PP and MIP and return a top list. This can create the appearance of having preferred one model over another when, in fact, the PP of each model is very similar. It is therefore important to check the distribution of the PP across the models. A low PP for the top model or very little variation across the top 5 or 10 models in the PP can be an indicator that there has been an error in the implementation of the algorithm or that the data being used is unsuitable for this type of estimation. In our application, the top two models, including glutamine and glucose and only glucose, respectively, had PP's of 0.377 and 0.311. In contrast, the third highest model had a PP of 0.051, showing much lower evidence of being the true causal model.

If the MVMR‐BMA does struggle to distinguish between different models, estimation could potentially be improved by restricting the set of exposures further. For example, by considering each sub‐class of metabolite separately. Although such an approach reduces the ‘hypothesis‐free’ aspect of the analysis, a ‘hypothesis‐free' approach is not possible when the genetic data is not sufficient to distinguish between noise and a true causal signal. The success of this approach will depend on the researchers' ability to reduce the exposure traits to a potentially meaningful set, where each set can be identified by the genetic instruments selected in the data.

Extensions to the MVMR‐BMA method that would also help to overcome this issue include incorporating prior probability to exposures, weak instrument robust estimation techniques, or by limiting the SNPs used for each model considered to those associated with the exposures included in that model. Neither approachhas been implemented within the algorithm so far, and each would have its own strengths and limitations. Therefore, further research is required to determine if these alterations would improve the results of the model when applied to data such as that considered here.

MVMR‐BMA reports a MIP, which is the probability that the final model includes that exposure. It can be hard to determine from a MIP whether an exposure should be selected as an exposure, as it is unclear what probability is high enough to result in inclusion of the exposure in the model. Empirical p‐values can be calculated for the MVMR‐BMA estimation. In connection with multiple testing procedures like false discovery rates, these empirical p‐values can be used to define which exposures should be selected as having a likely causal effect on the outcome (Levin et al. 2021; Zuber, Gill, et al. ²⁰²⁰).