Perform MR

library(TwoSampleMR)
library(ggplot2)

Introduction

Let’s continue with the example of BMI on CHD:

bmi_exp_dat <- extract_instruments(outcomes = 'ieu-a-2')
chd_out_dat <- extract_outcome_data(snps = bmi_exp_dat$SNP, outcomes = 'ieu-a-7')
dat <- harmonise_data(bmi_exp_dat, chd_out_dat)
#> Harmonising Body mass index || id:ieu-a-2 (ieu-a-2) and Coronary heart disease || id:ieu-a-7 (ieu-a-7)

Once the exposure and outcome data are harmonised, we have effects and standard errors for each instrument SNP available for the exposure and outcome traits. We can use this information to perform Mendelian randomisation. To do this, simply run:

res <- mr(dat)
#> Analysing 'ieu-a-2' on 'ieu-a-7'
res
#>   id.exposure id.outcome                              outcome
#> 1     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
#> 2     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
#> 3     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
#> 4     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
#> 5     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
#>                        exposure                    method nsnp         b
#> 1 Body mass index || id:ieu-a-2                  MR Egger   79 0.5024935
#> 2 Body mass index || id:ieu-a-2           Weighted median   79 0.3870065
#> 3 Body mass index || id:ieu-a-2 Inverse variance weighted   79 0.4459091
#> 4 Body mass index || id:ieu-a-2               Simple mode   79 0.3401554
#> 5 Body mass index || id:ieu-a-2             Weighted mode   79 0.3790910
#>           se         pval
#> 1 0.14396056 8.012590e-04
#> 2 0.07514313 2.601288e-07
#> 3 0.05898302 4.032020e-14
#> 4 0.16260431 3.969906e-02
#> 5 0.09975725 2.851076e-04

This returns a data frame of estimates of the causal effect of the exposure on the outcome for a range of different MR methods.

If there were multiple exposures against multiple outcomes in dat, the mr() function will perform each MR method for each combination of exposure-outcome traits.

MR methods

The list of available MR methods can be obtained:

mr_method_list()
#>                              obj
#> 1                  mr_wald_ratio
#> 2               mr_two_sample_ml
#> 3            mr_egger_regression
#> 4  mr_egger_regression_bootstrap
#> 5               mr_simple_median
#> 6             mr_weighted_median
#> 7   mr_penalised_weighted_median
#> 8                         mr_ivw
#> 9                  mr_ivw_radial
#> 10                    mr_ivw_mre
#> 11                     mr_ivw_fe
#> 12                mr_simple_mode
#> 13              mr_weighted_mode
#> 14         mr_weighted_mode_nome
#> 15           mr_simple_mode_nome
#> 16                       mr_raps
#> 17                       mr_sign
#> 18                        mr_uwr
#>                                                         name PubmedID
#> 1                                                 Wald ratio         
#> 2                                         Maximum likelihood         
#> 3                                                   MR Egger 26050253
#> 4                                       MR Egger (bootstrap) 26050253
#> 5                                              Simple median         
#> 6                                            Weighted median         
#> 7                                  Penalised weighted median         
#> 8                                  Inverse variance weighted         
#> 9                                                 IVW radial         
#> 10 Inverse variance weighted (multiplicative random effects)         
#> 11                 Inverse variance weighted (fixed effects)         
#> 12                                               Simple mode         
#> 13                                             Weighted mode         
#> 14                                      Weighted mode (NOME)         
#> 15                                        Simple mode (NOME)         
#> 16                      Robust adjusted profile score (RAPS)         
#> 17                                     Sign concordance test         
#> 18                                     Unweighted regression         
#>                                                    Description use_by_default
#> 1                                                                        TRUE
#> 2                                                                       FALSE
#> 3                                                                        TRUE
#> 4                                                                       FALSE
#> 5                                                                       FALSE
#> 6                                                                        TRUE
#> 7                                                                       FALSE
#> 8                                                                        TRUE
#> 9                                                                       FALSE
#> 10                                                                      FALSE
#> 11                                                                      FALSE
#> 12                                                                       TRUE
#> 13                                                                       TRUE
#> 14                                                                      FALSE
#> 15                                                                      FALSE
#> 16                                                                      FALSE
#> 17 Tests for concordance of signs between exposure and outcome          FALSE
#> 18                                     Doesn't use any weights          FALSE
#>    heterogeneity_test
#> 1               FALSE
#> 2                TRUE
#> 3                TRUE
#> 4               FALSE
#> 5               FALSE
#> 6               FALSE
#> 7               FALSE
#> 8                TRUE
#> 9                TRUE
#> 10              FALSE
#> 11              FALSE
#> 12              FALSE
#> 13              FALSE
#> 14              FALSE
#> 15              FALSE
#> 16              FALSE
#> 17              FALSE
#> 18               TRUE

To perform them, they can be specified in the mr() function, e.g. to only perform MR Egger regression and Inverse variance weighted methods,

mr(dat, method_list = c("mr_egger_regression", "mr_ivw"))
#> Analysing 'ieu-a-2' on 'ieu-a-7'
#>   id.exposure id.outcome                              outcome
#> 1     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
#> 2     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
#>                        exposure                    method nsnp         b
#> 1 Body mass index || id:ieu-a-2                  MR Egger   79 0.5024935
#> 2 Body mass index || id:ieu-a-2 Inverse variance weighted   79 0.4459091
#>           se        pval
#> 1 0.14396056 8.01259e-04
#> 2 0.05898302 4.03202e-14

By default, all the methods that are labelled TRUE in the use_by_default column are used by the mr() function.


Sensitivity analyses

Heterogeneity statistics

Some of the MR methods can also perform tests for heterogeneity. To obtain those statistics:

mr_heterogeneity(dat)
#>   id.exposure id.outcome                              outcome
#> 1     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
#> 2     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
#>                        exposure                    method        Q Q_df
#> 1 Body mass index || id:ieu-a-2                  MR Egger 143.3046   77
#> 2 Body mass index || id:ieu-a-2 Inverse variance weighted 143.6508   78
#>         Q_pval
#> 1 6.841585e-06
#> 2 8.728420e-06

As with the mr() function, the mr_heterogeneity() function can take an argument to only perform heterogeneity tests using specified methods, e.g.

mr_heterogeneity(dat, method_list = c("mr_egger_regression", "mr_ivw"))
#>   id.exposure id.outcome                              outcome
#> 1     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
#> 2     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
#>                        exposure                    method        Q Q_df
#> 1 Body mass index || id:ieu-a-2                  MR Egger 143.3046   77
#> 2 Body mass index || id:ieu-a-2 Inverse variance weighted 143.6508   78
#>         Q_pval
#> 1 6.841585e-06
#> 2 8.728420e-06

Horizontal pleiotropy

The intercept term in MR Egger regression can be a useful indication of whether directional horizontal pleiotropy is driving the results of an MR analysis. This can be obtained as follows:

mr_pleiotropy_test(dat)
#>   id.exposure id.outcome                              outcome
#> 1     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
#>                        exposure egger_intercept          se      pval
#> 1 Body mass index || id:ieu-a-2    -0.001719304 0.003985962 0.6674266

Single SNP analysis

To obtain the MR estimates using each of the SNPs singly we can do the following:

res_single <- mr_singlesnp(dat)

This returns a data.frame of results that is similar to the output from mr() except it performs the analysis multiple times for each exposure-outcome combination - each time using a different single SNP to perform the analysis.

The method used to perform the single SNP MR is the Wald ratio by default, though this can be changed, e.g. to use the fixed effects meta analysis method instead:

res_single <- mr_singlesnp(dat, single_method = "mr_meta_fixed")

The mr_singlesnp() function calculates the full MR using all available SNPs as well, and by default it uses the IVW and MR Egger methods. This can be specified as so:

res_single <- mr_singlesnp(dat, all_method = "mr_two_sample_ml")

will perform only the maximum likelihood method for the combined test.

Leave-one-out analysis

It is possible to perform a leave-one-out analysis, where the MR is performed again but leaving out each SNP in turn, to identify if a single SNP is driving the association.

res_loo <- mr_leaveoneout(dat)

By default the method used is the inverse variance weighted method, but this can be changed by using the method argument.


Plots

There are a few ways to visualise the results, listed below

Scatter plot

We can depict the relationship of the SNP effects on the exposure against the SNP effects on the outcome using a scatter plot.

res <- mr(dat)
#> Analysing 'ieu-a-2' on 'ieu-a-7'
p1 <- mr_scatter_plot(res, dat)

A scatter plot is created for each exposure-outcome test, and stored in p1 as a list of plots. For example, to plot the first scatter plot:

p1[[1]]

A scatter plot visualising the two-sample data points and the following fitted models; Inverse Variance Weighted, MR-Egger, Simple mode, Weighted median, and Weighted mode.

And to see how many plots there are:

length(p1)
#> [1] 1

Lines are drawn for each method used in mr(dat), the slope of the line corresponding to the estimated causal effect. To limit which lines are drawn, simply specify the desired methods, e.g. to only draw MR Egger and IVW:

res <- mr(dat, method_list = c("mr_egger_regression", "mr_ivw"))
#> Analysing 'ieu-a-2' on 'ieu-a-7'
p1 <- mr_scatter_plot(res, dat)

It is possible to save this plot using the ggsave() function from the ggplot2 package, e.g. to save as a pdf

ggsave(p1[[1]], file = "filename.pdf", width = 7, height = 7)

or save as a png

ggsave(p1[[1]], file = "filename.png", width = 7, height = 7)

See ?ggplot2::ggsave for more info.

Forest plot

Use the mr_forest_plot() function to compare the MR estimates using the different MR methods against the single SNP tests.

res_single <- mr_singlesnp(dat)
p2 <- mr_forest_plot(res_single)
p2[[1]]

A forest plot showing the estimated causal effects using each SNP separately, and the Inverse Variance Weighted and MR-Egger estimates using all the SNPs.

Here, the plot shows the causal effect as estimated using each of the SNPs on their own, and comparing against the causal effect as estimated using the methods that use all the SNPs.

To get plots that use different methods, specify them in the mr_singlesnp() function:

res_single <- mr_singlesnp(dat, all_method = c("mr_ivw", "mr_two_sample_ml"))
p2 <- mr_forest_plot(res_single)
p2[[1]]

An alternative forest plot showing the estimated causal effects using each SNP separately, and the Inverse Variance Weighted and Maximum Likelihood estimates using all the SNPs.

Leave-one-out plot

Use the mr_leaveoneout_plot() function to visualise the leave-one-out analysis:

res_loo <- mr_leaveoneout(dat)
p3 <- mr_leaveoneout_plot(res_loo)
p3[[1]]