Nothing
Mendelian randomization vignette"
In MendelianRandomization: Mendelian Randomization Package
MendelianRandomization package
MendelianRandomization is a package developed to carry out Mendelian randomization analyses using summarized genetic data in R. The package implements various methods to help assess whether a risk factor (also called an exposure) has a causal effect on an outcome.
Several additional references are available which describe the operation of this package, including:
-
Yavorska and Burgess (2017) "MendelianRandomization: an R package for performing Mendelian randomization analyses using summarized data", doi: 10.1093/ije/dyx034. This describes the basic operation of the package and core functions.
-
Broadbent et al (2020) "MendelianRandomization v0.5.0: updates to an R package for performing Mendelian randomization analyses using summarized data", doi: 10.12688/wellcomeopenres.16374.2. This describes updates to the package up to version 0.5.0.
-
Patel et al (2023) "MendelianRandomization v0.9.0: updates to an R package for performing Mendelian randomization analyses using summarized data". This describes updates to the package up to version 0.9.0.
library(MendelianRandomization) # loading the package
The Input
The package uses a special class called MRInput within the analyses in order to pass in all necessary information through one simple structure rather than inserting the data in parts. In order to make an MRInput object, one can do the following:
- specify values for each slot separately, or
- extract values from the PhenoScanner web-based database
We focus initially on the first option.
The MRInput object has the following "slots" :
- betaX and betaXse are both numeric vectors describing the associations of the genetic variants with the exposure. betaX are the beta-coefficients from univariable regression analyses of the exposure on each genetic variant in turn, and betaXse are the standard errors.
- betaY and betaYse are both numeric vectors describing the associations of the genetic variants with the outcome. betaY are the beta-coefficients from regression analyses of the outcome on each genetic variant in turn, and betaYse are the standard errors.
- correlation is a matrix with the signed correlations between the variants. If a correlation matrix is not provided, it is assumed that the variants are uncorrelated.
- exposure is a character string giving the name of the risk factor, e.g. LDL-cholesterol.
- outcome is a character string giving the name of the outcome, e.g. coronary heart disease. These inputs are only used in the graphing functions.
- snps is a character vector of the names of the various genetic variants (SNPs) in the dataset, e.g. rs12785878. It is not necessary to name the exposure, outcome, or SNPs, but these names are used in the graphing functions and may be helpful for keeping track of various analyses.
To generate the MRInput object slot by slot, one can use the mr_input() function :
MRInputObject <- mr_input(bx = ldlc,
bxse = ldlcse,
by = chdlodds,
byse = chdloddsse)
MRInputObject # example with uncorrelated variants
MRInputObject.cor <- mr_input(bx = calcium,
bxse = calciumse,
by = fastgluc,
byse = fastglucse,
corr = calc.rho)
MRInputObject.cor # example with correlated variants
It is not necessary for all the slots to be filled. For example, some of the methods do not require bxse to be specified; for example, the mr_ivw function will still run with bxse set to zeros. If the vectors bx, bxse, by, and byse are not of equal length, then an error will be reported. Note that the package does not implement any harmonization of associations to the same effect allele; this must be done by the user.
It is also possible to run the analysis using the syntax:
MRInputObject <- mr_input(ldlc, ldlcse, chdlodds, chdloddsse)
However, care must be taken in this case to give the vectors in the correct order (that is: bx, bxse, by, byse).
The data
Two sets of data are provided as part of this package:
- ldlc, ldlcse, hdlc, hdlcse, trig, trigse, chdlodds, chdloddsse: these are the associations (beta-coefficients and standard errors) of 28 genetic variants with LDL-cholesterol, HDL-cholesterol, triglycerides, and coronary heart disease (CHD) risk (associations with CHD risk are log odds ratios) taken from Waterworth et al (2011) "Genetic variants influencing circulating lipid levels and risk of coronary artery disease", doi: 10.1161/atvbaha.109.201020.
- calcium, calciumse, fastgluc, fastglucse: these are the associations (beta-coefficients and standard errors) of 7 genetic variants in the CASR gene region. These 7 variants are all correlated, and the correlation matrix is provided as calc.rho. These data were analysed in Burgess et al (2015) "Using published data in Mendelian randomization: a blueprint for efficient identification of causal risk factors", doi: 10.1007/s10654-015-0011-z.
Univariable estimation methods
The MendelianRandomization package supports several univariable estimation methods (that is, methods for a single exposure variable): the inverse-variance weighted method, the median-based method, the MR-Egger method, the mode-based method, the maximum likelihood method, the heterogeneity penalized method, the contamination mixture method, the lasso method, the debiased inverse-variance weighted method, the penalized inverse-variance weighted method, the constrained maximum likelihood method, and the principal components generalized method of moments (GMM) method. We describe these methods in turn.
Inverse-variance weighted method
The inverse-variance method is equivalent to the standard instrumental variable method using individual-level data, the two-stage least squares method. Either a fixed- or a random-effects analysis can be performed; the "default" option is a fixed-effect analysis when there are three variants or fewer, and a random-effects analysis otherwise. The robust option uses robust regression rather than standard regression in the analysis, and the penalized option downweights the contribution to the analysis of genetic variants with outlying (heterogeneous) variant-specific estimates. If a correlation matrix is provided in the MRInput object, then the correlated method is used by default (correl = TRUE), and the robust and penalized arguments are ignored.
The default options for constructing confidence intervals are based on a normal distribution and a 95% confidence level, however one can use the t-distribution (distribution = "t-dist") and alternative significance level if desired.
The default option for the weights is "simple", which corresponds to an inverse-variance weighted meta-analysis of the ratio estimates by/bx using first-order standard errors byse/bx. Alternatively, weights can be set to "delta", to include second-order terms of the delta expansion. This makes use of psi, which is the correlation for each variant between its association with the exposure and with the outcome.
The mr_ivw function automatically calculates an approximation to the first-stage F statistic from regression of the exposure on the genetic variants. This is an important measure of instrument strength. When the correlation matrix is specified, calculation of the estimated F statistic accounts for the correlation between variants. When a correlation matrix is specified, a warning message appears to remind the user of the importance of correct harmonization of the correlation matrix to the same effect and other alleles as the genetic associations with the exposure and outcome.
IVWObject <- mr_ivw(MRInputObject,
model = "default",
robust = FALSE,
penalized = FALSE,
correl = FALSE,
weights = "simple",
psi = 0,
distribution = "normal",
alpha = 0.05)
IVWObject <- mr_ivw(mr_input(bx = ldlc, bxse = ldlcse,
by = chdlodds, byse = chdloddsse))
IVWObject
IVWObject.correl <- mr_ivw(MRInputObject.cor,
model = "default",
correl = TRUE,
distribution = "normal",
alpha = 0.05)
IVWObject.correl <- mr_ivw(mr_input(bx = calcium, bxse = calciumse,
by = fastgluc, byse = fastglucse, corr = calc.rho))
IVWObject.correl
Median-based method
The median-based method calculates a median of the variant-specific estimates from the ratio method for each genetic variant individually. The default option is to calculate a weighted median using the inverse-variance weights. Alternatively, one can calculate a simple (unweighted) median, or a weighted median using penalization of weights for heterogeneous variants. Since the calculation of standard error requires bootstrapping, the number of bootstrap iterations can be varied. The random seed is set automatically so that results are reproducible; however, the value of the seed can be changed if required. Once the function has been completed, the random seed is reset to its previous value (this is true for all functions in the package).
The median-based method requires data on at least 3 genetic variants. Variants must be uncorrelated.
WeightedMedianObject <- mr_median(MRInputObject,
weighting = "weighted",
distribution = "normal",
alpha = 0.05,
iterations = 10000,
seed = 314159265)
WeightedMedianObject <- mr_median(mr_input(bx = ldlc, bxse = ldlcse,
by = chdlodds, byse = chdloddsse))
WeightedMedianObject
SimpleMedianObject <- mr_median(mr_input(bx = ldlc, bxse = ldlcse,
by = chdlodds, byse = chdloddsse), weighting = "simple")
SimpleMedianObject
MR-Egger method
The MR-Egger method is implemented here using a random-effects model only. The robust and penalized options are the same as for the inverse-variance weighted method. The method can be used for both correlated and uncorrelated sets of variants. Confidence intervals can be constructed either using a normal distribution (distribution = "normal", the default option), or a t-distribution (distribution = "t-dist").
With a t-distribution, in case of under-dispersion (the estimated residual standard error in the regression model is less than 1), confidence intervals and p-values use either a t-distribution with no correction for under-dispersion, or a normal distribution with the residual standard error set to 1 -- whichever is wider. This means that under-dispersion is not doubly penalized by setting the residual standard error to 1 and using a t-distribution, but also that the confidence intervals are not narrower (p-values not more extreme) than those using a fixed-effect model.
The MR-Egger method requires data on at least 3 genetic variants. Variants are permitted to be correlated.
EggerObject <- mr_egger(MRInputObject,
robust = FALSE,
penalized = FALSE,
correl = FALSE,
distribution = "normal",
alpha = 0.05)
EggerObject <- mr_egger(mr_input(bx = ldlc, bxse = ldlcse,
by = chdlodds, byse = chdloddsse))
EggerObject
EggerObject.corr <- mr_egger(MRInputObject.cor,
correl = TRUE,
distribution = "normal",
alpha = 0.05)
EggerObject.corr <- mr_egger(mr_input(bx = calcium, bxse = calciumse,
by = fastgluc, byse = fastglucse, corr = calc.rho))
EggerObject.corr
Maximum likelihood method
An alternative estimation method is the maximum likelihood method, introduced in Burgess et al (2013) "Mendelian randomization analysis with multiple genetic variants using summarized data", doi: 10.1002/gepi.21758.
The method has two main advantages over the IVW method: it allows for uncertainty in the genetic associations with the exposure (which is ignored in the IVW method using simple weights), and it allows for genetic associations with the exposure and with the outcome for each variant to be correlated. This correlation arises if the samples for the associations with the exposure and the outcome overlap. In a strict two-sample Mendelian randomization setting, the correlation parameter psi is set to zero (the default option). If the associations are estimated in the same individuals (complete sample overlap), then the observed correlation between the exposure and the outcome is a reasonable proposal for the value of psi.
MaxLikObject <- mr_maxlik(MRInputObject,
model = "default",
correl = FALSE,
psi = 0,
distribution = "normal",
alpha = 0.05)
MaxLikObject <- mr_maxlik(mr_input(bx = ldlc, bxse = ldlcse,
by = chdlodds, byse = chdloddsse))
MaxLikObject
MaxLikObject.corr <- mr_maxlik(mr_input(bx = calcium, bxse = calciumse,
by = fastgluc, byse = fastglucse, corr = calc.rho))
MaxLikObject.corr
Mode-based estimation method
An additional method for robust estimation that gives consistent when a plurality of genetic variants is valid is the mode-based estimation method, introduced in Hartwig, Davey Smith and Bowden (2017) "Robust inference in summary data Mendelian randomization via the zero modal pleiotropy assumption", doi: 10.1093/ije/dyx102.
While the median-based method calculates the variant-specific estimates and then obtains the median of these estimates, and the inverse-variance weighted estimate is a weighted mean of the variant-specific estimates, the mode-based method obtains the "mode" of these estimates. However, in finite samples, no two estimates will be exactly the same, so a mode does not exist. The method proceeds by constructing a kernel-weighted density of the variant-specific estimates, and taking the maximum point of this density as the point estimate. A confidence interval is obtained by bootstrapping.
Several options can be specified: whether the mode should be "weighted" or "unweighted" (default is weighted), whether the standard errors should be calculated using the "simple" or "delta" formula (default is delta), the value of the bandwidth parameter multiplication factor (default is 1), the random seed, and number of iterations in the bootstrap procedure.
MBEObject <- mr_mbe(MRInputObject,
weighting = "weighted",
stderror = "delta",
phi = 1,
seed = 314159265,
iterations = 10000,
distribution = "normal",
alpha = 0.05)
MBEObject <- mr_mbe(mr_input(bx = ldlc, bxse = ldlcse,
by = chdlodds, byse = chdloddsse))
MBEObject
Heterogeneity-penalized method
Another method for robust estimation that also gives consistent when a plurality of genetic variants is valid is the heterogeneity-penalized model-averaging method, introduced in Burgess et al (2018) "Modal-based estimation via heterogeneity-penalized weighting: model averaging for consistent and efficient estimation in Mendelian randomization when a plurality of candidate instruments are valid", doi: 10.1093/ije/dyy08010.1101/175372.
The heterogeneity-penalized method uses the same consistency criterion as the mode-based estimation method, but evaluates the modal estimate in a different way. It proceeds by evaluating weights for all subsets of genetic variants (excluding the null set and singletons). Subsets receive greater weight if they include more variants, but are severely downweighted if the variants in the subset have heterogeneous variant-specific estimates. As such, the method will identify the subset with the largest number (by weight) of variants having similar variant-specific estimates.
Confidence intervals are evaluated by calculating a log-likelihood function, and finding all points within a given vertical distance of the maximum of the log-likelihood function (which is the overall estimate). As such, if the log-likelihood function is multimodal, then the confidence interval may include multiple disjoint ranges. This may indicate the presence of multiple causal mechanisms by which the exposure may influence the outcome with different magnitudes of causal effect. As the confidence interval is determined by a grid search, care must be taken when chosing the minimum (CIMin) and maximum (CIMax) values in the search, as well as the step size (CIStep). The default values will not be suitable for all applications.
HetPenObject <- mr_hetpen(MRInputObject,
prior = 0.5,
CIMin = -1,
CIMax = 1,
CIStep = 0.001,
alpha = 0.05)
As the method evaluates estimates and weights for each subset of variants, the method is quite computationally expensive to run, as the complexity doubles with each additional variant. We run the method for a subset of 10 genetic variants:
HetPenObject <- mr_hetpen(mr_input(bx = ldlc[1:10], bxse = ldlcse[1:10],
by = chdlodds[1:10], byse = chdloddsse[1:10]), CIMin = -1, CIMax = 5, CIStep = 0.01)
HetPenObject
As an example of a multimodal confidence interval:
bcrp =c(0.160, 0.236, 0.149, 0.09, 0.079, 0.072, 0.047, 0.069)
bcrpse =c(0.006, 0.009, 0.006, 0.005, 0.005, 0.005, 0.006, 0.011)
bchd =c(0.0237903, -0.1121942, -0.0711906, -0.030848, 0.0479207, 0.0238895,
0.005528, 0.0214852)
bchdse =c(0.0149064, 0.0303084, 0.0150552, 0.0148339, 0.0143077, 0.0145478,
0.0160765, 0.0255237)
HetPenObject.multimode <- mr_hetpen(mr_input(bx = bcrp, bxse = bcrpse,
by = bchd, byse = bchdse))
HetPenObject.multimode
Other univariable estimation methods
Several other estimation methods are available, which we do not describe in detail, as their operation is similar to methods already described above. We would encourage interested analysts to read the help files for each method, which describe the inputs and outputs of each method.
The contamination mixture method is introduced in Burgess et al (2020) "A robust and efficient method for Mendelian randomization with hundreds of genetic variants", doi: 10.1038/s41467-019-14156-4. It is similar in spirit to the heterogeneity-penalized method, but much more computationally efficient.
ConMixObject <- mr_conmix(MRInputObject,
psi = 0,
CIMin = NA,
CIMax = NA,
CIStep = 0.01,
alpha = 0.05)
ConMixObject <- mr_conmix(mr_input(bx = ldlc, bxse = ldlcse,
by = chdlodds, byse = chdloddsse))
ConMixObject
The lasso method is introduced in Rees et al (2019) "Robust methods in Mendelian randomization via penalization of heterogeneous causal estimates", doi: 10.1371/journal.pone.0222362. It is an outlier-robust method, which excludes variants with heterogeneous variant-specific estimates using lasso penalization.
LassoObject <- mr_lasso(MRInputObject,
distribution = "normal",
alpha = 0.05,
lambda = numeric(0))
LassoObject <- mr_lasso(mr_input(bx = ldlc, bxse = ldlcse,
by = chdlodds, byse = chdloddsse))
LassoObject
The debiased inverse-variance weighted method is a variation on the inverse-variance weighted method that is more robust to bias due to weak instruments and winner's curse. It is described in Ting, Shao, Kang (2021) "Debiased inverse-variance weighted estimator in two-sample summary-data Mendelian randomization", doi: 10.1214/20-aos2027.
DIVWObject <- mr_divw(MRInputObject,
over.dispersion = TRUE,
alpha = 0.05,
diagnostics = FALSE)
DIVWObject <- mr_divw(mr_input(bx = ldlc, bxse = ldlcse,
by = chdlodds, byse = chdloddsse))
DIVWObject
The penalized inverse-variance weighted method is a further extension of the inverse-variance weighted and debiased inverse-variance weighted methods to deal with selection of variants more explicitly. It is described in Wang et al (2022) "A novel penalized inverse-variance weighted estimator for Mendelian randomization with applications to COVID-19 outcomes", doi: 10.1111/biom.13732.
PIVWObject <- mr_pivw(MRInputObject,
over.dispersion = TRUE,
delta = 0,
sel.pval = NULL,
Boot.Fieller = TRUE,
alpha = 0.05)
PIVWObject <- mr_pivw(mr_input(bx = ldlc, bxse = ldlcse,
by = chdlodds, byse = chdloddsse))
PIVWObject
The constrained maximum likelihood method is a robust method that can account for violation of any of the three instrumental variable assumptions under mild assumptions. In a maximum likelihood framework, this method constrains the number of invalid instruments with horizontal pleiotropy. The number of invalid instruments is asymptotically consistently selected by the Bayesian information criterion. It is described in Xue, Shen, Pan (2021) "Constrained maximum likelihood-based Mendelian randomization robust to both correlated and uncorrelated pleiotropic effects", doi: 10.1016/j.ajhg.2021.05.014.
Note that this method requires the sample size used to calculate the genetic associations with the exposure and/or outcome to be specified (if these differ, then it is recommended to use the lower value; see reference for details).
cMLObject <- mr_cML(MRInputObject,
MA = TRUE,
DP = TRUE,
K_vec = 0:(length(object@betaX)-2),
random_start = 0,
num_pert = 200,
random_start_pert = 0,
maxit = 100,
random_seed = 314,
n,
Alpha = 0.05)
cMLObject <- mr_cML(mr_input(bx = ldlc, bxse = ldlcse,
by = chdlodds, byse = chdloddsse), n = 17723)
cMLObject
The principal components GMM method is an implementation of the GMM method with summarized data that is designed for use when performing Mendelian randomization using genetic variants from a single gene region. As an alternative to pruning and clumping approaches, which take a large number of variants from a gene region (potentially hundreds or thousands) and select a small number of uncorrelated (or weakly correlated) variants, the principal components approach performs dimension reduction on the full set of variant associations. The principal components are then used as instrument variables, rather than the individual genetic variants. The method is described in Patel et al (2023) "Robust use of phenotypic heterogeneity at drug target genes for mechanistic insights: application of cis-multivariable Mendelian randomization to GLP1R gene region", doi: 10.1101/2023.07.20.23292958.
Note that this method requires the sample size used to calculate the genetic associations with the exposure (nx) and outcome (ny) to be specified. The number of principal components includes as instruments can either be set directly, by setting r, or indirectly, by setting thres. The default option thres = 0.999 includes enough principal components to explain 99.9% of the variance in the weighted variant correlation matrix.
pcGMMObject <- mr_pcgmm(MRInputObject.cor,
nx,
ny,
r = NULL,
thres = 0.999,
robust = TRUE,
alpha = 0.05)
pcGMMObject <- mr_pcgmm(mr_input(bx = calcium, bxse = calciumse,
by = fastgluc, byse = fastglucse, corr = calc.rho),
nx=6351, ny=133010)
pcGMMObject
Multivariable Mendelian randomization
An analytic approach when genetic variants are associated with multiple risk factors is multivariable Mendelian randomization, introduced in Burgess and Thompson (2015) "Multivariable Mendelian randomization: the use of pleiotropic genetic variants to estimate causal effects", doi: 10.1093/aje/kwu283.
There are two contexts in which we envision the method being used. First, when there is a set of related risk factors, such as lipid fractions. It is difficult to find genetic predictors of HDL-cholesterol that do not also predict LDL-cholesterol and/or triglycerides. Multivariable Mendelian randomization allows genetic variants to be associated with all the risk factors in the model provided that they do not influence the outcome via other pathways. Secondly, when there is a network of risk factors, typically a primary risk factor and a secondary risk factor or mediator. In both cases, multivariable estimates reflect direct causal effects - the result of intervening on the risk factor under analysis while keeping other risk factors in the model fixed.
As the multivariable method requires genetic associations with multiple risk factors, a different input function is needed: mr_mvinput, which creates an MRMVInput object:
MVMRInputObject <- mr_mvinput(bx = cbind(ldlc, hdlc, trig),
bxse = cbind(ldlcse, hdlcse, trigse),
by = chdlodds,
byse = chdloddsse)
MVMRInputObject
MVMRInputObject.cor <- mr_mvinput(bx = cbind(ldlc, hdlc, trig),
bxse = cbind(ldlcse, hdlcse, trigse),
by = chdlodds,
byse = chdloddsse,
correlation = diag(length(ldlc)))
The mr_mvivw command performs an extension of the inverse-variance weighted method for multivariable Mendelian randomization, implementing multivariable weighted linear regression (the standard IVW method method performs univariable weighted linear regression) with the intercept term set to zero. As with the standard IVW method, a correlation matrix can be specified. Multivariable Mendelian randomization requires the number of genetic variants to exceed the number of risk factors; if this is not the case, the method will return an error.
MVIVWObject <- mr_mvivw(MVMRInputObject,
model = "default",
correl = FALSE,
correl.x = NULL,
nx = NA,
distribution = "normal",
alpha = 0.05)
MVIVWObject <- mr_mvivw(MVMRInputObject)
MVIVWObject
Conditional F statistics: The nx option specifies the sample size(s) for the genetic associations with the exposures (either a single value if these are all equal, or a vector if they differ). The correl.x option specifies the correlation between the genetic associations with the exposures. If the genetic assocations with the exposures are estimates in separate datasets, then this correlation matrix should be set to the identity matrix (this is the default option). This matrix cannot be calculated from summarized data; if it is unknown, then a sensitivity analysis for its value is recommended. The nx and correl.x options are not used in the calculation of estimates from the multivariable inverse-variance weighted method method, but only in the calculation of conditional F statistics. Conditional F statistics are the relevant measure of instrument strength for multivariable Mendelian randomization. For more detail, see Sanderson, Spiller, Bowden (2021) "Testing and correcting for weak and pleiotropic instruments in two-sample multivariable Mendelian randomization", doi: 10.1002/sim.9133.
MVIVWObject.condF <- mr_mvivw(MVMRInputObject, nx = 17723)
MVIVWObject.condF
Other univariable estimation methods
In addition to the multivariable inverse-variance weighted method, several other multivariable estimation methods are available, which generally are multivariable extensions of univariable methods: multivariable MR-Egger, multivariable median-based method, multivariable lasso method, multivariable constrained maximum likelihood, multivariable GMM, and multivariable principal components GMM.
MVEggerObject <- mr_mvegger(MVMRInputObject)
MVEggerObject
MVMedianObject <- mr_mvmedian(MVMRInputObject)
MVMedianObject
MVLassoObject <- mr_mvlasso(MVMRInputObject)
MVLassoObject
MVcMLObject <- mr_mvcML(MVMRInputObject, n = 17723)
MVcMLObject
MVGMMObject <- mr_mvgmm(MVMRInputObject, nx=rep(17723,3), ny=17723)
MVGMMObject
MVpcGMMObject <- mr_mvpcgmm(MVMRInputObject.cor, nx=rep(17723,3), ny=17723)
MVpcGMMObject
Note that the mr_mvpcgmm example does not use variants from a single gene region, and so is provided to demonstrate that the code works, rather than to illustrate a recommended use case.
Summaries of multiple methods
The mr_allmethods function is provided to easily compare results (Estimate, Standard Error, 95% CI, and p-value) from multiple methods. One can look at results from a wide range of methods (method = "all"), or a limited set of results by setting method to "egger", "ivw", or "median". The final option is "main", which gives results from the simple median, weighted median, IVW, and MR-Egger methods only.
MRInputObject <- mr_input(bx = ldlc,
bxse = ldlcse,
by = chdlodds,
byse = chdloddsse)
MRAllObject_all <- mr_allmethods(MRInputObject, method = "all")
MRAllObject_all
MRAllObject_egger <- mr_allmethods(MRInputObject, method = "egger")
MRAllObject_egger
MRAllObject_main <- mr_allmethods(MRInputObject, method = "main")
MRAllObject_main
Graphical summaries of results
Applied to MRInput object - display of data
The mr_plot function has two different functionalities. First, if the function is applied to an MRInput object, then the output is an interactive graph that can be used to explore the associations of the different genetic variants and look for outliers, which may represent pleiotropic variants.
The syntax is:
mr_plot(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse),
error = TRUE, orientate = FALSE, line = "ivw")
An interactive graph does not reproduce well in a vignette, so we encourage readers to input this code for themselves. The interactive graph allows the user to pinpoint outliers easily; when the user mouses over one of the points, the name of the variant is shown.
The option error = TRUE plots error bars (95% confidence intervals) for the associations with the exposure and with the outcome. The option orientate = TRUE sets all the associations with the exposure to be positive, and re-orientates the associations with the outcome if needed. This option is encouraged for the MR-Egger method (as otherwise points having negative associations with the exposure can appear to be far from the regression line), although by default it is set to FALSE. The line option can be set to "ivw" (to show the inverse-variance weighted estimate) or to "egger" (to show the MR-Egger estimate).
A static version of this graph is also available, by setting interactive = FALSE:
mr_plot(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse),
error = TRUE, orientate = FALSE, line = "ivw", interactive = FALSE)
This version of the graph is less useful for detecting outliers, but easier to save as a file and share with colleagues.
Outliers can be identified using the labels = TRUE option:
mr_plot(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse),
error = TRUE, orientate = FALSE, line = "ivw", interactive = FALSE, labels = TRUE)
The resulting graph is quite ugly, but it is easy to identify the individual points.
Applied to MRMVInput object - display of data
The mr_plot command can also be applied to an MRMVInput object created using the mr_mvinput() function. In this case, the horizontal axis does not present the associations with any single exposure, but the fitted value from the multivariable IVW method, representing the expected genetic association with the outcome based on the genetic associations with the exposures.
mr_plot(MVMRInputObject, interactive = FALSE)
Applied to MRAll object - comparison of estimates
Finally, if the mr_plot function is applied to the output of the mr_allmethods function, estimates from the different methods can be compared graphically.
mr_plot(MRAllObject_all)
mr_plot(MRAllObject_egger)
mr_plot(mr_allmethods(mr_input(bx = hdlc, bxse = hdlcse,
by = chdlodds, byse = chdloddsse)))
We see that estimates from all methods are similar when LDL-cholesterol is the risk factor, but the MR-Egger estimates differ substantially when HDL-cholesterol is the risk factor.
Other graphical functions
Other graphical functions are the mr_forest command, which can either provide a forest plot of the variant-specific estimates:
mr_forest(MRInputObject, ordered=TRUE)
Or the comparison of various estimates from different methods:
mr_forest(MRInputObject,
methods = c("ivw", "median", "wmedian", "egger", "maxlik", "mbe", "conmix"),
snp_estimates = FALSE)
The mr_funnel command, which provides a funnel plot of the variant-specific estimates against their precision (that is, the reciprocal of the standard error):
mr_funnel(mr_input(bx = ldlc[1:8], bxse = ldlcse[1:8],
by = chdlodds[1:8], byse = chdloddsse[1:8]))
(Note we only plot the first 8 variants here)
And the mr_loo command, which provides the leave-one-out estimates - that is, the estimate based on all variants except variant 1, then the based on all variants except variant 2, and so on.
mr_loo(MRInputObject)
These plots are all ggplot objects, and so can be edited using commands from the ggplot2 package. For example, to set the x-axis to run from -5 to +5:
library(ggplot2)
forest = mr_forest(mr_input(ldlc, ldlcse, chdlodds, chdloddsse))
forest2 = forest + coord_cartesian(xlim=c(-5,5))
forest2
Extracting association estimates from PhenoScanner
The PhenoScanner bioinformatic tool (http://phenoscanner.medschl.cam.ac.uk) is a curated database of publicly available results from large-scale genetic association studies. The database currently contains over 65 billion associations and association results and over 150 million unique genetic variants, mostly single nucleotide polymorphisms.
PhenoScanner can be called directly from the MendelianRandomization package using the pheno_input() function. This creates an MRInput function, which can be directly used as an input to any of the estimation functions. For example:
mr_ivw(pheno_input(snps=c("rs12916", "rs2479409",
"rs217434", "rs1367117",
"rs4299376", "rs629301",
"rs4420638", "rs6511720"),
exposure = "Low density lipoprotein",
pmidE = "24097068",
ancestryE = "European",
outcome = "Coronary artery disease",
pmidO = "26343387",
ancestryO = "Mixed"))
(We do not implement this code here, as it requires a connection to the internet, and hence produces an error if an internet connection cannot be found. But please copy it and try it for yourself!)
In order to obtain the relevant summary estimates, run the pheno_input() function with:
- snps is a character vector giving the rsid identifiers of the genetic variants.
- exposure is a character vector giving the name of the risk factor.
- pmidE is the PubMed ID of the paper where the association estimates with the exposure were first published.
- ancestryE is the ancestry of the participants on whom the association estimates with the exposure were estimated. (For some traits and PubMed IDs, results are given for multiple ancestries.) Usually, ancestry is "European" or "Mixed".
- outcome is a character vector giving the name of the outcome.
- pmidO is the PubMed ID of the paper where the association estimates with the outcome were first published.
- ancestryO is the ancestry of the participants on whom the association estimates with the exposure were estimated.
We note that the spelling of the exposure and outcome, the PubMed ID, and the ancestry information need to correspond exactly to the values in the PhenoScanner dataset. If these are not spelled exactly as in the PhenoScanner dataset (including upper/lower case), the association estimates will not be found.
Final note of caution
Particularly with the development of Mendelian randomization with summarized data, two-sample Mendelian randomization (where associations with the risk factor and with the outcome are taken from separate datasets), and bioinformatic tools for obtaining and analysing summarized data (including this package), Mendelian randomization is becoming increasingly accessible as a tool for answering epidemiological questions. This is undoubtably a Good Thing.
However, it is important to remember that the difficult part of a Mendelian randomization analysis is not the computational method, but deciding what should go into the analysis: which risk factor, which outcome, and (most importantly) which genetic variants. Hopefully, the availability of these tools will enable less attention to be paid to the mechanics of the analysis, and more attention to these choices.
The note of caution is that tools that make Mendelian randomization simple to perform run the risk of encouraging large numbers of speculative analyses to be performed in an unprincipled way. It is important that Mendelian randomization is not performed in a way that avoids critical thought. In releasing this package, the hope is that it will lead to more comprehensive and more reproducible causal inferences from Mendelian randomization, and not simply add more noise to the literature.
Try the MendelianRandomization package in your browser
Any scripts or data that you put into this service are public.
MendelianRandomization documentation built on May 29, 2024, 11:36 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
MendelianRandomization package
MendelianRandomization is a package developed to carry out Mendelian randomization analyses using summarized genetic data in R. The package implements various methods to help assess whether a risk factor (also called an exposure) has a causal effect on an outcome.
Several additional references are available which describe the operation of this package, including:
-
Yavorska and Burgess (2017) "MendelianRandomization: an R package for performing Mendelian randomization analyses using summarized data", doi: 10.1093/ije/dyx034. This describes the basic operation of the package and core functions.
-
Broadbent et al (2020) "MendelianRandomization v0.5.0: updates to an R package for performing Mendelian randomization analyses using summarized data", doi: 10.12688/wellcomeopenres.16374.2. This describes updates to the package up to version 0.5.0.
-
Patel et al (2023) "MendelianRandomization v0.9.0: updates to an R package for performing Mendelian randomization analyses using summarized data". This describes updates to the package up to version 0.9.0.
library(MendelianRandomization) # loading the package
The Input
The package uses a special class called MRInput within the analyses in order to pass in all necessary information through one simple structure rather than inserting the data in parts. In order to make an MRInput object, one can do the following:
- specify values for each slot separately, or
- extract values from the PhenoScanner web-based database
We focus initially on the first option.
The MRInput object has the following "slots" :
- betaX and betaXse are both numeric vectors describing the associations of the genetic variants with the exposure. betaX are the beta-coefficients from univariable regression analyses of the exposure on each genetic variant in turn, and betaXse are the standard errors.
- betaY and betaYse are both numeric vectors describing the associations of the genetic variants with the outcome. betaY are the beta-coefficients from regression analyses of the outcome on each genetic variant in turn, and betaYse are the standard errors.
- correlation is a matrix with the signed correlations between the variants. If a correlation matrix is not provided, it is assumed that the variants are uncorrelated.
- exposure is a character string giving the name of the risk factor, e.g. LDL-cholesterol.
- outcome is a character string giving the name of the outcome, e.g. coronary heart disease. These inputs are only used in the graphing functions.
- snps is a character vector of the names of the various genetic variants (SNPs) in the dataset, e.g. rs12785878. It is not necessary to name the exposure, outcome, or SNPs, but these names are used in the graphing functions and may be helpful for keeping track of various analyses.
To generate the MRInput object slot by slot, one can use the mr_input() function :
MRInputObject <- mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse) MRInputObject # example with uncorrelated variants MRInputObject.cor <- mr_input(bx = calcium, bxse = calciumse, by = fastgluc, byse = fastglucse, corr = calc.rho) MRInputObject.cor # example with correlated variants
It is not necessary for all the slots to be filled. For example, some of the methods do not require bxse to be specified; for example, the mr_ivw function will still run with bxse set to zeros. If the vectors bx, bxse, by, and byse are not of equal length, then an error will be reported. Note that the package does not implement any harmonization of associations to the same effect allele; this must be done by the user.
It is also possible to run the analysis using the syntax:
MRInputObject <- mr_input(ldlc, ldlcse, chdlodds, chdloddsse)
However, care must be taken in this case to give the vectors in the correct order (that is: bx, bxse, by, byse).
The data
Two sets of data are provided as part of this package:
- ldlc, ldlcse, hdlc, hdlcse, trig, trigse, chdlodds, chdloddsse: these are the associations (beta-coefficients and standard errors) of 28 genetic variants with LDL-cholesterol, HDL-cholesterol, triglycerides, and coronary heart disease (CHD) risk (associations with CHD risk are log odds ratios) taken from Waterworth et al (2011) "Genetic variants influencing circulating lipid levels and risk of coronary artery disease", doi: 10.1161/atvbaha.109.201020.
- calcium, calciumse, fastgluc, fastglucse: these are the associations (beta-coefficients and standard errors) of 7 genetic variants in the CASR gene region. These 7 variants are all correlated, and the correlation matrix is provided as calc.rho. These data were analysed in Burgess et al (2015) "Using published data in Mendelian randomization: a blueprint for efficient identification of causal risk factors", doi: 10.1007/s10654-015-0011-z.
Univariable estimation methods
The MendelianRandomization package supports several univariable estimation methods (that is, methods for a single exposure variable): the inverse-variance weighted method, the median-based method, the MR-Egger method, the mode-based method, the maximum likelihood method, the heterogeneity penalized method, the contamination mixture method, the lasso method, the debiased inverse-variance weighted method, the penalized inverse-variance weighted method, the constrained maximum likelihood method, and the principal components generalized method of moments (GMM) method. We describe these methods in turn.
Inverse-variance weighted method
The inverse-variance method is equivalent to the standard instrumental variable method using individual-level data, the two-stage least squares method. Either a fixed- or a random-effects analysis can be performed; the "default" option is a fixed-effect analysis when there are three variants or fewer, and a random-effects analysis otherwise. The robust option uses robust regression rather than standard regression in the analysis, and the penalized option downweights the contribution to the analysis of genetic variants with outlying (heterogeneous) variant-specific estimates. If a correlation matrix is provided in the MRInput object, then the correlated method is used by default (correl = TRUE), and the robust and penalized arguments are ignored.
The default options for constructing confidence intervals are based on a normal distribution and a 95% confidence level, however one can use the t-distribution (distribution = "t-dist") and alternative significance level if desired.
The default option for the weights is "simple", which corresponds to an inverse-variance weighted meta-analysis of the ratio estimates by/bx using first-order standard errors byse/bx. Alternatively, weights can be set to "delta", to include second-order terms of the delta expansion. This makes use of psi, which is the correlation for each variant between its association with the exposure and with the outcome.
The mr_ivw function automatically calculates an approximation to the first-stage F statistic from regression of the exposure on the genetic variants. This is an important measure of instrument strength. When the correlation matrix is specified, calculation of the estimated F statistic accounts for the correlation between variants. When a correlation matrix is specified, a warning message appears to remind the user of the importance of correct harmonization of the correlation matrix to the same effect and other alleles as the genetic associations with the exposure and outcome.
IVWObject <- mr_ivw(MRInputObject, model = "default", robust = FALSE, penalized = FALSE, correl = FALSE, weights = "simple", psi = 0, distribution = "normal", alpha = 0.05) IVWObject <- mr_ivw(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse)) IVWObject IVWObject.correl <- mr_ivw(MRInputObject.cor, model = "default", correl = TRUE, distribution = "normal", alpha = 0.05) IVWObject.correl <- mr_ivw(mr_input(bx = calcium, bxse = calciumse, by = fastgluc, byse = fastglucse, corr = calc.rho)) IVWObject.correl
Median-based method
The median-based method calculates a median of the variant-specific estimates from the ratio method for each genetic variant individually. The default option is to calculate a weighted median using the inverse-variance weights. Alternatively, one can calculate a simple (unweighted) median, or a weighted median using penalization of weights for heterogeneous variants. Since the calculation of standard error requires bootstrapping, the number of bootstrap iterations can be varied. The random seed is set automatically so that results are reproducible; however, the value of the seed can be changed if required. Once the function has been completed, the random seed is reset to its previous value (this is true for all functions in the package).
The median-based method requires data on at least 3 genetic variants. Variants must be uncorrelated.
WeightedMedianObject <- mr_median(MRInputObject, weighting = "weighted", distribution = "normal", alpha = 0.05, iterations = 10000, seed = 314159265) WeightedMedianObject <- mr_median(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse)) WeightedMedianObject SimpleMedianObject <- mr_median(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse), weighting = "simple") SimpleMedianObject
MR-Egger method
The MR-Egger method is implemented here using a random-effects model only. The robust and penalized options are the same as for the inverse-variance weighted method. The method can be used for both correlated and uncorrelated sets of variants. Confidence intervals can be constructed either using a normal distribution (distribution = "normal", the default option), or a t-distribution (distribution = "t-dist").
With a t-distribution, in case of under-dispersion (the estimated residual standard error in the regression model is less than 1), confidence intervals and p-values use either a t-distribution with no correction for under-dispersion, or a normal distribution with the residual standard error set to 1 -- whichever is wider. This means that under-dispersion is not doubly penalized by setting the residual standard error to 1 and using a t-distribution, but also that the confidence intervals are not narrower (p-values not more extreme) than those using a fixed-effect model.
The MR-Egger method requires data on at least 3 genetic variants. Variants are permitted to be correlated.
EggerObject <- mr_egger(MRInputObject, robust = FALSE, penalized = FALSE, correl = FALSE, distribution = "normal", alpha = 0.05) EggerObject <- mr_egger(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse)) EggerObject EggerObject.corr <- mr_egger(MRInputObject.cor, correl = TRUE, distribution = "normal", alpha = 0.05) EggerObject.corr <- mr_egger(mr_input(bx = calcium, bxse = calciumse, by = fastgluc, byse = fastglucse, corr = calc.rho)) EggerObject.corr
Maximum likelihood method
An alternative estimation method is the maximum likelihood method, introduced in Burgess et al (2013) "Mendelian randomization analysis with multiple genetic variants using summarized data", doi: 10.1002/gepi.21758.
The method has two main advantages over the IVW method: it allows for uncertainty in the genetic associations with the exposure (which is ignored in the IVW method using simple weights), and it allows for genetic associations with the exposure and with the outcome for each variant to be correlated. This correlation arises if the samples for the associations with the exposure and the outcome overlap. In a strict two-sample Mendelian randomization setting, the correlation parameter psi is set to zero (the default option). If the associations are estimated in the same individuals (complete sample overlap), then the observed correlation between the exposure and the outcome is a reasonable proposal for the value of psi.
MaxLikObject <- mr_maxlik(MRInputObject, model = "default", correl = FALSE, psi = 0, distribution = "normal", alpha = 0.05) MaxLikObject <- mr_maxlik(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse)) MaxLikObject MaxLikObject.corr <- mr_maxlik(mr_input(bx = calcium, bxse = calciumse, by = fastgluc, byse = fastglucse, corr = calc.rho)) MaxLikObject.corr
Mode-based estimation method
An additional method for robust estimation that gives consistent when a plurality of genetic variants is valid is the mode-based estimation method, introduced in Hartwig, Davey Smith and Bowden (2017) "Robust inference in summary data Mendelian randomization via the zero modal pleiotropy assumption", doi: 10.1093/ije/dyx102.
While the median-based method calculates the variant-specific estimates and then obtains the median of these estimates, and the inverse-variance weighted estimate is a weighted mean of the variant-specific estimates, the mode-based method obtains the "mode" of these estimates. However, in finite samples, no two estimates will be exactly the same, so a mode does not exist. The method proceeds by constructing a kernel-weighted density of the variant-specific estimates, and taking the maximum point of this density as the point estimate. A confidence interval is obtained by bootstrapping.
Several options can be specified: whether the mode should be "weighted" or "unweighted" (default is weighted), whether the standard errors should be calculated using the "simple" or "delta" formula (default is delta), the value of the bandwidth parameter multiplication factor (default is 1), the random seed, and number of iterations in the bootstrap procedure.
MBEObject <- mr_mbe(MRInputObject, weighting = "weighted", stderror = "delta", phi = 1, seed = 314159265, iterations = 10000, distribution = "normal", alpha = 0.05) MBEObject <- mr_mbe(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse)) MBEObject
Heterogeneity-penalized method
Another method for robust estimation that also gives consistent when a plurality of genetic variants is valid is the heterogeneity-penalized model-averaging method, introduced in Burgess et al (2018) "Modal-based estimation via heterogeneity-penalized weighting: model averaging for consistent and efficient estimation in Mendelian randomization when a plurality of candidate instruments are valid", doi: 10.1093/ije/dyy08010.1101/175372.
The heterogeneity-penalized method uses the same consistency criterion as the mode-based estimation method, but evaluates the modal estimate in a different way. It proceeds by evaluating weights for all subsets of genetic variants (excluding the null set and singletons). Subsets receive greater weight if they include more variants, but are severely downweighted if the variants in the subset have heterogeneous variant-specific estimates. As such, the method will identify the subset with the largest number (by weight) of variants having similar variant-specific estimates.
Confidence intervals are evaluated by calculating a log-likelihood function, and finding all points within a given vertical distance of the maximum of the log-likelihood function (which is the overall estimate). As such, if the log-likelihood function is multimodal, then the confidence interval may include multiple disjoint ranges. This may indicate the presence of multiple causal mechanisms by which the exposure may influence the outcome with different magnitudes of causal effect. As the confidence interval is determined by a grid search, care must be taken when chosing the minimum (CIMin) and maximum (CIMax) values in the search, as well as the step size (CIStep). The default values will not be suitable for all applications.
HetPenObject <- mr_hetpen(MRInputObject, prior = 0.5, CIMin = -1, CIMax = 1, CIStep = 0.001, alpha = 0.05)
As the method evaluates estimates and weights for each subset of variants, the method is quite computationally expensive to run, as the complexity doubles with each additional variant. We run the method for a subset of 10 genetic variants:
HetPenObject <- mr_hetpen(mr_input(bx = ldlc[1:10], bxse = ldlcse[1:10], by = chdlodds[1:10], byse = chdloddsse[1:10]), CIMin = -1, CIMax = 5, CIStep = 0.01) HetPenObject
As an example of a multimodal confidence interval:
bcrp =c(0.160, 0.236, 0.149, 0.09, 0.079, 0.072, 0.047, 0.069) bcrpse =c(0.006, 0.009, 0.006, 0.005, 0.005, 0.005, 0.006, 0.011) bchd =c(0.0237903, -0.1121942, -0.0711906, -0.030848, 0.0479207, 0.0238895, 0.005528, 0.0214852) bchdse =c(0.0149064, 0.0303084, 0.0150552, 0.0148339, 0.0143077, 0.0145478, 0.0160765, 0.0255237) HetPenObject.multimode <- mr_hetpen(mr_input(bx = bcrp, bxse = bcrpse, by = bchd, byse = bchdse)) HetPenObject.multimode
Other univariable estimation methods
Several other estimation methods are available, which we do not describe in detail, as their operation is similar to methods already described above. We would encourage interested analysts to read the help files for each method, which describe the inputs and outputs of each method.
The contamination mixture method is introduced in Burgess et al (2020) "A robust and efficient method for Mendelian randomization with hundreds of genetic variants", doi: 10.1038/s41467-019-14156-4. It is similar in spirit to the heterogeneity-penalized method, but much more computationally efficient.
ConMixObject <- mr_conmix(MRInputObject, psi = 0, CIMin = NA, CIMax = NA, CIStep = 0.01, alpha = 0.05) ConMixObject <- mr_conmix(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse)) ConMixObject
The lasso method is introduced in Rees et al (2019) "Robust methods in Mendelian randomization via penalization of heterogeneous causal estimates", doi: 10.1371/journal.pone.0222362. It is an outlier-robust method, which excludes variants with heterogeneous variant-specific estimates using lasso penalization.
LassoObject <- mr_lasso(MRInputObject, distribution = "normal", alpha = 0.05, lambda = numeric(0)) LassoObject <- mr_lasso(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse)) LassoObject
The debiased inverse-variance weighted method is a variation on the inverse-variance weighted method that is more robust to bias due to weak instruments and winner's curse. It is described in Ting, Shao, Kang (2021) "Debiased inverse-variance weighted estimator in two-sample summary-data Mendelian randomization", doi: 10.1214/20-aos2027.
DIVWObject <- mr_divw(MRInputObject, over.dispersion = TRUE, alpha = 0.05, diagnostics = FALSE) DIVWObject <- mr_divw(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse)) DIVWObject
The penalized inverse-variance weighted method is a further extension of the inverse-variance weighted and debiased inverse-variance weighted methods to deal with selection of variants more explicitly. It is described in Wang et al (2022) "A novel penalized inverse-variance weighted estimator for Mendelian randomization with applications to COVID-19 outcomes", doi: 10.1111/biom.13732.
PIVWObject <- mr_pivw(MRInputObject, over.dispersion = TRUE, delta = 0, sel.pval = NULL, Boot.Fieller = TRUE, alpha = 0.05) PIVWObject <- mr_pivw(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse)) PIVWObject
The constrained maximum likelihood method is a robust method that can account for violation of any of the three instrumental variable assumptions under mild assumptions. In a maximum likelihood framework, this method constrains the number of invalid instruments with horizontal pleiotropy. The number of invalid instruments is asymptotically consistently selected by the Bayesian information criterion. It is described in Xue, Shen, Pan (2021) "Constrained maximum likelihood-based Mendelian randomization robust to both correlated and uncorrelated pleiotropic effects", doi: 10.1016/j.ajhg.2021.05.014.
Note that this method requires the sample size used to calculate the genetic associations with the exposure and/or outcome to be specified (if these differ, then it is recommended to use the lower value; see reference for details).
cMLObject <- mr_cML(MRInputObject, MA = TRUE, DP = TRUE, K_vec = 0:(length(object@betaX)-2), random_start = 0, num_pert = 200, random_start_pert = 0, maxit = 100, random_seed = 314, n, Alpha = 0.05)
cMLObject <- mr_cML(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse), n = 17723) cMLObject
The principal components GMM method is an implementation of the GMM method with summarized data that is designed for use when performing Mendelian randomization using genetic variants from a single gene region. As an alternative to pruning and clumping approaches, which take a large number of variants from a gene region (potentially hundreds or thousands) and select a small number of uncorrelated (or weakly correlated) variants, the principal components approach performs dimension reduction on the full set of variant associations. The principal components are then used as instrument variables, rather than the individual genetic variants. The method is described in Patel et al (2023) "Robust use of phenotypic heterogeneity at drug target genes for mechanistic insights: application of cis-multivariable Mendelian randomization to GLP1R gene region", doi: 10.1101/2023.07.20.23292958.
Note that this method requires the sample size used to calculate the genetic associations with the exposure (nx) and outcome (ny) to be specified. The number of principal components includes as instruments can either be set directly, by setting r, or indirectly, by setting thres. The default option thres = 0.999 includes enough principal components to explain 99.9% of the variance in the weighted variant correlation matrix.
pcGMMObject <- mr_pcgmm(MRInputObject.cor, nx, ny, r = NULL, thres = 0.999, robust = TRUE, alpha = 0.05)
pcGMMObject <- mr_pcgmm(mr_input(bx = calcium, bxse = calciumse, by = fastgluc, byse = fastglucse, corr = calc.rho), nx=6351, ny=133010) pcGMMObject
Multivariable Mendelian randomization
An analytic approach when genetic variants are associated with multiple risk factors is multivariable Mendelian randomization, introduced in Burgess and Thompson (2015) "Multivariable Mendelian randomization: the use of pleiotropic genetic variants to estimate causal effects", doi: 10.1093/aje/kwu283.
There are two contexts in which we envision the method being used. First, when there is a set of related risk factors, such as lipid fractions. It is difficult to find genetic predictors of HDL-cholesterol that do not also predict LDL-cholesterol and/or triglycerides. Multivariable Mendelian randomization allows genetic variants to be associated with all the risk factors in the model provided that they do not influence the outcome via other pathways. Secondly, when there is a network of risk factors, typically a primary risk factor and a secondary risk factor or mediator. In both cases, multivariable estimates reflect direct causal effects - the result of intervening on the risk factor under analysis while keeping other risk factors in the model fixed.
As the multivariable method requires genetic associations with multiple risk factors, a different input function is needed: mr_mvinput, which creates an MRMVInput object:
MVMRInputObject <- mr_mvinput(bx = cbind(ldlc, hdlc, trig), bxse = cbind(ldlcse, hdlcse, trigse), by = chdlodds, byse = chdloddsse) MVMRInputObject MVMRInputObject.cor <- mr_mvinput(bx = cbind(ldlc, hdlc, trig), bxse = cbind(ldlcse, hdlcse, trigse), by = chdlodds, byse = chdloddsse, correlation = diag(length(ldlc)))
The mr_mvivw command performs an extension of the inverse-variance weighted method for multivariable Mendelian randomization, implementing multivariable weighted linear regression (the standard IVW method method performs univariable weighted linear regression) with the intercept term set to zero. As with the standard IVW method, a correlation matrix can be specified. Multivariable Mendelian randomization requires the number of genetic variants to exceed the number of risk factors; if this is not the case, the method will return an error.
MVIVWObject <- mr_mvivw(MVMRInputObject, model = "default", correl = FALSE, correl.x = NULL, nx = NA, distribution = "normal", alpha = 0.05) MVIVWObject <- mr_mvivw(MVMRInputObject) MVIVWObject
Conditional F statistics: The nx option specifies the sample size(s) for the genetic associations with the exposures (either a single value if these are all equal, or a vector if they differ). The correl.x option specifies the correlation between the genetic associations with the exposures. If the genetic assocations with the exposures are estimates in separate datasets, then this correlation matrix should be set to the identity matrix (this is the default option). This matrix cannot be calculated from summarized data; if it is unknown, then a sensitivity analysis for its value is recommended. The nx and correl.x options are not used in the calculation of estimates from the multivariable inverse-variance weighted method method, but only in the calculation of conditional F statistics. Conditional F statistics are the relevant measure of instrument strength for multivariable Mendelian randomization. For more detail, see Sanderson, Spiller, Bowden (2021) "Testing and correcting for weak and pleiotropic instruments in two-sample multivariable Mendelian randomization", doi: 10.1002/sim.9133.
MVIVWObject.condF <- mr_mvivw(MVMRInputObject, nx = 17723) MVIVWObject.condF
Other univariable estimation methods
In addition to the multivariable inverse-variance weighted method, several other multivariable estimation methods are available, which generally are multivariable extensions of univariable methods: multivariable MR-Egger, multivariable median-based method, multivariable lasso method, multivariable constrained maximum likelihood, multivariable GMM, and multivariable principal components GMM.
MVEggerObject <- mr_mvegger(MVMRInputObject) MVEggerObject MVMedianObject <- mr_mvmedian(MVMRInputObject) MVMedianObject MVLassoObject <- mr_mvlasso(MVMRInputObject) MVLassoObject MVcMLObject <- mr_mvcML(MVMRInputObject, n = 17723) MVcMLObject MVGMMObject <- mr_mvgmm(MVMRInputObject, nx=rep(17723,3), ny=17723) MVGMMObject MVpcGMMObject <- mr_mvpcgmm(MVMRInputObject.cor, nx=rep(17723,3), ny=17723) MVpcGMMObject
Note that the mr_mvpcgmm example does not use variants from a single gene region, and so is provided to demonstrate that the code works, rather than to illustrate a recommended use case.
Summaries of multiple methods
The mr_allmethods function is provided to easily compare results (Estimate, Standard Error, 95% CI, and p-value) from multiple methods. One can look at results from a wide range of methods (method = "all"), or a limited set of results by setting method to "egger", "ivw", or "median". The final option is "main", which gives results from the simple median, weighted median, IVW, and MR-Egger methods only.
MRInputObject <- mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse) MRAllObject_all <- mr_allmethods(MRInputObject, method = "all") MRAllObject_all MRAllObject_egger <- mr_allmethods(MRInputObject, method = "egger") MRAllObject_egger MRAllObject_main <- mr_allmethods(MRInputObject, method = "main") MRAllObject_main
Graphical summaries of results
Applied to MRInput object - display of data
The mr_plot function has two different functionalities. First, if the function is applied to an MRInput object, then the output is an interactive graph that can be used to explore the associations of the different genetic variants and look for outliers, which may represent pleiotropic variants.
The syntax is:
mr_plot(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse), error = TRUE, orientate = FALSE, line = "ivw")
An interactive graph does not reproduce well in a vignette, so we encourage readers to input this code for themselves. The interactive graph allows the user to pinpoint outliers easily; when the user mouses over one of the points, the name of the variant is shown.
The option error = TRUE plots error bars (95% confidence intervals) for the associations with the exposure and with the outcome. The option orientate = TRUE sets all the associations with the exposure to be positive, and re-orientates the associations with the outcome if needed. This option is encouraged for the MR-Egger method (as otherwise points having negative associations with the exposure can appear to be far from the regression line), although by default it is set to FALSE. The line option can be set to "ivw" (to show the inverse-variance weighted estimate) or to "egger" (to show the MR-Egger estimate).
A static version of this graph is also available, by setting interactive = FALSE:
mr_plot(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse), error = TRUE, orientate = FALSE, line = "ivw", interactive = FALSE)
This version of the graph is less useful for detecting outliers, but easier to save as a file and share with colleagues.
Outliers can be identified using the labels = TRUE option:
mr_plot(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse), error = TRUE, orientate = FALSE, line = "ivw", interactive = FALSE, labels = TRUE)
The resulting graph is quite ugly, but it is easy to identify the individual points.
Applied to MRMVInput object - display of data
The mr_plot command can also be applied to an MRMVInput object created using the mr_mvinput() function. In this case, the horizontal axis does not present the associations with any single exposure, but the fitted value from the multivariable IVW method, representing the expected genetic association with the outcome based on the genetic associations with the exposures.
mr_plot(MVMRInputObject, interactive = FALSE)
Applied to MRAll object - comparison of estimates
Finally, if the mr_plot function is applied to the output of the mr_allmethods function, estimates from the different methods can be compared graphically.
mr_plot(MRAllObject_all)
mr_plot(MRAllObject_egger)
mr_plot(mr_allmethods(mr_input(bx = hdlc, bxse = hdlcse, by = chdlodds, byse = chdloddsse)))
We see that estimates from all methods are similar when LDL-cholesterol is the risk factor, but the MR-Egger estimates differ substantially when HDL-cholesterol is the risk factor.
Other graphical functions
Other graphical functions are the mr_forest command, which can either provide a forest plot of the variant-specific estimates:
mr_forest(MRInputObject, ordered=TRUE)
Or the comparison of various estimates from different methods:
mr_forest(MRInputObject, methods = c("ivw", "median", "wmedian", "egger", "maxlik", "mbe", "conmix"), snp_estimates = FALSE)
The mr_funnel command, which provides a funnel plot of the variant-specific estimates against their precision (that is, the reciprocal of the standard error):
mr_funnel(mr_input(bx = ldlc[1:8], bxse = ldlcse[1:8], by = chdlodds[1:8], byse = chdloddsse[1:8]))
(Note we only plot the first 8 variants here)
And the mr_loo command, which provides the leave-one-out estimates - that is, the estimate based on all variants except variant 1, then the based on all variants except variant 2, and so on.
mr_loo(MRInputObject)
These plots are all ggplot objects, and so can be edited using commands from the ggplot2 package. For example, to set the x-axis to run from -5 to +5:
library(ggplot2) forest = mr_forest(mr_input(ldlc, ldlcse, chdlodds, chdloddsse)) forest2 = forest + coord_cartesian(xlim=c(-5,5)) forest2
Extracting association estimates from PhenoScanner
The PhenoScanner bioinformatic tool (http://phenoscanner.medschl.cam.ac.uk) is a curated database of publicly available results from large-scale genetic association studies. The database currently contains over 65 billion associations and association results and over 150 million unique genetic variants, mostly single nucleotide polymorphisms.
PhenoScanner can be called directly from the MendelianRandomization package using the pheno_input() function. This creates an MRInput function, which can be directly used as an input to any of the estimation functions. For example:
mr_ivw(pheno_input(snps=c("rs12916", "rs2479409", "rs217434", "rs1367117", "rs4299376", "rs629301", "rs4420638", "rs6511720"), exposure = "Low density lipoprotein", pmidE = "24097068", ancestryE = "European", outcome = "Coronary artery disease", pmidO = "26343387", ancestryO = "Mixed"))
(We do not implement this code here, as it requires a connection to the internet, and hence produces an error if an internet connection cannot be found. But please copy it and try it for yourself!)
In order to obtain the relevant summary estimates, run the pheno_input() function with: - snps is a character vector giving the rsid identifiers of the genetic variants. - exposure is a character vector giving the name of the risk factor. - pmidE is the PubMed ID of the paper where the association estimates with the exposure were first published. - ancestryE is the ancestry of the participants on whom the association estimates with the exposure were estimated. (For some traits and PubMed IDs, results are given for multiple ancestries.) Usually, ancestry is "European" or "Mixed". - outcome is a character vector giving the name of the outcome. - pmidO is the PubMed ID of the paper where the association estimates with the outcome were first published. - ancestryO is the ancestry of the participants on whom the association estimates with the exposure were estimated.
We note that the spelling of the exposure and outcome, the PubMed ID, and the ancestry information need to correspond exactly to the values in the PhenoScanner dataset. If these are not spelled exactly as in the PhenoScanner dataset (including upper/lower case), the association estimates will not be found.
Final note of caution
Particularly with the development of Mendelian randomization with summarized data, two-sample Mendelian randomization (where associations with the risk factor and with the outcome are taken from separate datasets), and bioinformatic tools for obtaining and analysing summarized data (including this package), Mendelian randomization is becoming increasingly accessible as a tool for answering epidemiological questions. This is undoubtably a Good Thing.
However, it is important to remember that the difficult part of a Mendelian randomization analysis is not the computational method, but deciding what should go into the analysis: which risk factor, which outcome, and (most importantly) which genetic variants. Hopefully, the availability of these tools will enable less attention to be paid to the mechanics of the analysis, and more attention to these choices.
The note of caution is that tools that make Mendelian randomization simple to perform run the risk of encouraging large numbers of speculative analyses to be performed in an unprincipled way. It is important that Mendelian randomization is not performed in a way that avoids critical thought. In releasing this package, the hope is that it will lead to more comprehensive and more reproducible causal inferences from Mendelian randomization, and not simply add more noise to the literature.
Try the MendelianRandomization package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.