Skip to contents

Perform Mendelian Randomization Bayesian Model Averaging

Source: R/mr_bma.R
mr_bma.Rd

This function takes harmonized exposure/outcome data and runs Mendelian Randomization Bayesian Model Averaging to identify and prioritize relationships between the exposures and outcome.

Usage

mr_bma(
  harmonized_data,
  prior_prob = 0.5,
  prior_sigma = 0.5,
  top = 10,
  kmin,
  kmax,
  remove_outliers = TRUE,
  remove_influential = TRUE,
  calculate_p = FALSE,
  nrepeat = 1e+05
)

Arguments

harmonized_data

(list) Harmonized data of exposures/outcomes of interest, generated using TwoSampleMR::mv_harmonise_data()

prior_prob

(numeric) Prior probability

prior_sigma

(numeric) Variance of the prior probability

top

(numeric) Number of top models to return in results

kmin

(numeric) Minimum model size. By default the model will consider all combinations of exposures (eg. kmin = kmax), but may become computationally infeasible with >12 exposures. If kmin < kmax, then a stochastic search is performed.

kmax

(numeric) Maximum model size. By default the model will consider all combinations of exposures (eg. kmin = kmax), but may become computationally infeasible with >12 exposures. If kmin < kmax, then a stochastic search is performed.

remove_outliers

(logical) Remove outlier variants based on Q-statistic

remove_influential

(logical) Remove influential variants based on Cook's distance

calculate_p

(logical) Use empirical permutation with Nyholt procedure of effective tests to estimate p-values

nrepeat

(numeric) Number of permutations used when estimating Nyholt-corrected False Discovery Rate

Value

A list containing the results of MR-BMA analyses. The list includes:

  • model_best = A tibble containing a list of the top models

  • mip_table = A tibble containing the marginal inclusion probabilities of each risk factor

  • mrbma_output = Raw output from the summarymvMR_SSS function

  • influential_res = Diagnostic plots representing the detection of influential variants

  • outlier_res = Diagnostic plot representing the detection of outlier variants

Examples

# Extract genetic instruments
lipid_exposures <- TwoSampleMR::mv_extract_exposures(id_exposure = c("ieu-a-299", "ieu-a-300", "ieu-a-302"))
#> API: public: http://gwas-api.mrcieu.ac.uk/
#> Please look at vignettes for options on running this locally if you need to run many instances of this command.
#> Clumping 1, 214 variants, using EUR population reference
#> Removing 70 of 214 variants due to LD with other variants or absence from LD reference panel
#> Extracting data for 144 SNP(s) from 3 GWAS(s)
#> Harmonising HDL cholesterol || id:ieu-a-299 (ieu-a-299) and LDL cholesterol || id:ieu-a-300 (ieu-a-300)
#> Harmonising HDL cholesterol || id:ieu-a-299 (ieu-a-299) and Triglycerides || id:ieu-a-302 (ieu-a-302)

# Extract corresponding outcome data
cad_outcome <- TwoSampleMR::extract_outcome_data(snps = lipid_exposures$SNP, outcomes = "ebi-a-GCST005195")
#> Extracting data for 144 SNP(s) from 1 GWAS(s)
#> Finding proxies for 2 SNPs in outcome ebi-a-GCST005195
#> Extracting data for 2 SNP(s) from 1 GWAS(s)

# Generate harmonized dataset
lipids_cad_harmonized <- TwoSampleMR::mv_harmonise_data(lipid_exposures, cad_outcome)
#> Harmonising HDL cholesterol || id:ieu-a-299 (ieu-a-299) and Coronary artery disease || id:ebi-a-GCST005195 (ebi-a-GCST005195)
#> Removing the following SNPs for being palindromic with intermediate allele frequencies:
#> rs7534572

# Run MR-BMA
mr_bma_res <- mr_bma(lipids_cad_harmonized, calculate_p = TRUE, nrepeat = 1000)
#> ℹ Preparing input
#> ✔ Preparing input [87ms]
#> 
#> ℹ Running MR-BMA
#> ✔ Running MR-BMA [50ms]
#> 
#> ℹ Removing influential variants
#> ✔ Removing influential variants [98ms]
#> 
#> ℹ Removing outliers
#> ✔ Removing outliers [106ms]
#> 
#> ℹ Estimating empirical p-value using 1000 permutations
#> ✔ Estimating empirical p-value using 1000 permutations [2.9s]
#> 
#> ℹ Performing Nyholt correction for effective number of tests
#> ✔ Performing Nyholt correction for effective number of tests [32ms]
#> 

# Output best models
mr_bma_res$model_best
#> # A tibble: 7 × 3
#>   rf_combination                posterior_prob causal_estimate  
#>   <chr>                                  <dbl> <chr>            
#> 1 ieu-a-300,ieu-a-302                    0.843 0.365,0.219      
#> 2 ieu-a-299,ieu-a-300,ieu-a-302          0.156 -0.063,0.361,0.19
#> 3 ieu-a-299,ieu-a-300                    0     -0.131,0.376     
#> 4 ieu-a-300                              0     0.393            
#> 5 ieu-a-302                              0     0.311            
#> 6 ieu-a-299,ieu-a-302                    0     -0.105,0.26      
#> 7 ieu-a-299                              0     -0.202           

# Output marginal inclusion probabilities for each risk factor
mr_bma_res$mip_table
#> # A tibble: 3 × 6
#>   exposure  mace                 mip               p_val           p_fdr p_nyh…¹
#>   <chr>     <chr>                <chr>             <chr>           <chr> <chr>  
#> 1 ieu-a-299 -0.00987190503160796 0.15651065273458  0.927072927072… 0.92… 1      
#> 2 ieu-a-300 0.364151410155237    1                 0.000999000999… 0.00… 0.0019…
#> 3 ieu-a-302 0.214705609094214    0.999726011149142 0.000999000999… 0.00… 0.0019…
#> # … with abbreviated variable name ¹​p_nyholt