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.
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.
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
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:
To obtain the MR estimates using each of the SNPs singly we can do the following:
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:
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:
will perform only the maximum likelihood method for the combined test.
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.
By default the method used is the inverse variance weighted method,
but this can be changed by using the method
argument.
There are a few ways to visualise the results, listed below
We can depict the relationship of the SNP effects on the exposure against the SNP effects on the outcome using a scatter plot.
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:
And to see how many plots there are:
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
or save as a png
See ?ggplot2::ggsave
for more info.
Use the mr_forest_plot()
function to compare the MR
estimates using the different MR methods against the single SNP
tests.
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]]