Univariable estimation methods
Source:vignettes/Univariable_MR_Methods.Rmd
Univariable_MR_Methods.RmdThe 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.
IVW method: Example with uncorrelated variants
MRInputObject <- mr_input(bx = ldlc,
bxse = ldlcse,
by = chdlodds,
byse = chdloddsse)
MRInputObject # example with uncorrelated variants## SNP exposure.beta exposure.se outcome.beta outcome.se
## 1 snp_1 0.026 0.004 0.0677 0.0286
## 2 snp_2 -0.044 0.004 -0.1625 0.0300
## 3 snp_3 -0.038 0.004 -0.1054 0.0310
## 4 snp_4 -0.023 0.003 -0.0619 0.0243
## 5 snp_5 -0.017 0.003 -0.0834 0.0222
## 6 snp_6 -0.031 0.006 -0.1278 0.0667
## [ reached 'max' / getOption("max.print") -- omitted 22 rows ]
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##
## Inverse-variance weighted method
## (variants uncorrelated, random-effect model)
##
## Number of Variants : 28
##
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value
## IVW 2.834 0.530 1.796, 3.873 0.000
## ------------------------------------------------------------------
## Residual standard error = 1.920
## Heterogeneity test statistic (Cochran's Q) = 99.5304 on 27 degrees of freedom, (p-value = 0.0000). I^2 = 72.9%.
## F statistic = 28.0.
MRInputObject.cor <- mr_input(bx = calcium,
bxse = calciumse,
by = fastgluc,
byse = fastglucse,
corr = calc.rho)IVW method: example with correlated variants
MRInputObject.cor # example with correlated variants## SNP exposure.beta exposure.se outcome.beta outcome.se
## 1 snp_1 0.00625 0.00233 0.02805 0.0122
## 2 snp_2 0.00590 0.00338 0.00953 0.0198
## 3 snp_3 0.01822 0.00318 0.03646 0.0173
## 4 snp_4 0.00598 0.00233 0.01049 0.0119
## 5 snp_5 0.00825 0.00229 0.02357 0.0122
## 6 snp_6 0.00651 0.00352 0.00204 0.0179
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##
## Inverse-variance weighted method
## (variants correlated, random-effect model)
##
## Number of Variants : 6
##
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value
## IVW 2.245 0.643 0.984, 3.505 0.000
## ------------------------------------------------------------------
## Residual standard error = 0.641
## Residual standard error is set to 1 in calculation of confidence interval when its estimate is less than 1.
## Heterogeneity test statistic (Cochran's Q) = 2.0530 on 5 degrees of freedom, (p-value = 0.8418). I^2 = 0.0%.
## F statistic = 11.5.
##
## (Estimates with correlated variants are sensitive to the signs in the correlation matrix
## - please ensure that your correlations are expressed with respect to the same effect alleles as your summarized association estimates.)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 ##
## Weighted median method
##
## Number of Variants : 28
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value
## Weighted median method 2.683 0.419 1.862, 3.504 0.000
## ------------------------------------------------------------------
SimpleMedianObject <- mr_median(mr_input(bx = ldlc, bxse = ldlcse,
by = chdlodds, byse = chdloddsse), weighting = "simple")
SimpleMedianObject##
## Simple median method
##
## Number of Variants : 28
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value
## Simple median method 1.755 0.740 0.305, 3.205 0.018
## ------------------------------------------------------------------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##
## MR-Egger method
## (variants uncorrelated, random-effect model)
##
## Number of Variants = 28
##
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value
## MR-Egger 3.253 0.770 1.743, 4.762 0.000
## (intercept) -0.011 0.015 -0.041, 0.018 0.451
## ------------------------------------------------------------------
## Residual Standard Error : 1.935
## Heterogeneity test statistic = 97.3975 on 26 degrees of freedom, (p-value = 0.0000)
## I^2_GX statistic: 91.9%
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##
## MR-Egger method
## (variants correlated, random-effect model)
##
## Number of Variants = 6
##
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value
## MR-Egger 1.302 1.518 -1.672, 4.276 0.391
## (intercept) 0.009 0.013 -0.017, 0.035 0.493
## ------------------------------------------------------------------
## Residual Standard Error : 0.629
## Residual standard error is set to 1 in calculation of confidence interval when its estimate is less than 1.
## Heterogeneity test statistic = 1.5825 on 4 degrees of freedom, (p-value = 0.8119)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##
## Maximum-likelihood method
## (variants uncorrelated, random-effect model)
##
## Number of Variants : 28
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value
## MaxLik 3.225 0.569 2.110, 4.340 0.000
## ------------------------------------------------------------------
## Residual standard error = 1.793
## Heterogeneity test statistic = 86.8145 on 27 degrees of freedom, (p-value = 0.0000)
MaxLikObject.corr <- mr_maxlik(mr_input(bx = calcium, bxse = calciumse,
by = fastgluc, byse = fastglucse, corr = calc.rho))
MaxLikObject.corr##
## Maximum-likelihood method
## (variants correlated, random-effect model)
##
## Number of Variants : 6
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value
## MaxLik 2.303 0.709 0.914, 3.692 0.001
## ------------------------------------------------------------------
## Residual standard error = 0.588
## Residual standard error is set to 1 in calculation of confidence interval when its estimate is less than 1.
## Heterogeneity test statistic = 1.7277 on 5 degrees of freedom, (p-value = 0.8854)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##
## Mode-based method of Hartwig et al
## (weighted, delta standard errors [not assuming NOME], bandwidth factor = 1)
##
## Number of Variants : 28
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value
## MBE 2.946 1.005 0.977, 4.916 0.003
## ------------------------------------------------------------------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##
## Heterogeneity-penalized method
## (Prior probability of instrument validity = 0.5)
##
## Number of Variants : 10
## ------------------------------------------------------------------
## Method Estimate 95% CI
## HetPen 2.85 1.87, 3.90
## ------------------------------------------------------------------
## Note: confidence interval is a single range of values.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##
## Heterogeneity-penalized method
## (Prior probability of instrument validity = 0.5)
##
## Number of Variants : 8
## ------------------------------------------------------------------
## Method Estimate 95% CI
## HetPen -0.437 -0.599, -0.256
## 0.014, 0.449
## ------------------------------------------------------------------
## Note: confidence interval contains multiple ranges of values.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.
Contamination mixture 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##
## Contamination mixture method
## (Standard deviation of invalid estimands = 66.17905)
##
## Number of Variants : 28
## ------------------------------------------------------------------
## Method Estimate 95% CI p-value
## ConMix 2.73 2.00, 3.43 4.79e-15
## ------------------------------------------------------------------
## Note: confidence interval is a single range of values.Lasso method
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##
## MR-Lasso method
##
## Number of variants : 28
## Number of valid instruments : 26
## Tuning parameter : 0.7757552
## ------------------------------------------------------------------
## Exposure Estimate Std Error 95% CI p-value
## exposure 2.671 0.431 1.827, 3.515 0.000
## ------------------------------------------------------------------Debiased inverse-variance weighted method
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##
## Debiased inverse-variance weighted method
## (Over.dispersion:TRUE)
##
## Number of Variants : 28
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value Condition
## dIVW 2.940 0.531 1.900, 3.980 0.000 142.920
## ------------------------------------------------------------------Penalized inverse-variance weighted method
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##
## Penalized inverse-variance weighted method
##
## Over dispersion: TRUE
## Bootstrapping Fieller: TRUE
## Penalty parameter (lambda): 1
## IV selection threshold (delta): 0
## Number of variants : 28
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value Condition
## pIVW 2.932 0.528 1.922, 4.001 0.000 142.920
## ------------------------------------------------------------------Constrained maximum likelihood method
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##
## Constrained maximum likelihood method (MRcML)
## Number of Variants: 28
## Results for: cML-MA-BIC-DP
## ------------------------------------------------------------------
## Method Estimate SE Pvalue 95% CI
## cML-MA-BIC-DP 2.821 0.388 0.000 [2.060,3.582]
## ------------------------------------------------------------------Principal components GMM method
The principal components GMM method is an implementation of the generalized method of moments (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##
## Univariable principal components generalized method of moments (PC-GMM) method
##
## Number of principal components used : 4
##
## Robust model with overdispersion heterogeneity.
##
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value F-stat
## PC-GMM 2.476 0.837 0.836, 4.116 0.003 12.7
## ------------------------------------------------------------------
##
## Overdispersion heterogeneity parameter estimate = -10.49837
##
## Heterogeneity test statistic = 1.6842Conditional likelihood ratio (CLR) method
The mr_clr function calculates confidence intervals based on inverting the conditional likelihood ratio and other identification-robust tests.
In weak instrument settings, usual inference based on point estimates and standard errors may not be accurate. This method calculates confidence intervals based on inverting identification-robust tests proposed in Wang and Kang (“Weak-instrument robust tests in two-sample summary-data Mendelian randomization”, S. Wang and H. Kang, Biometrics, 2021.) that provide valid inferences regardless of instrument strength.
This includes conditional likelihood ratio (CLR), Kleibergen (K), and Anderson and Rubin (AR) tests.
Evidence from the econometrics literature suggests that CLR inference is the best option in terms of power under a wide range of settings.
Please note that these methods do not provide point estimates, only confidence intervals. While most examples provide a confidence interval that is a single range of values, in some cases the confidence interval may comprise multiple ranges of values. In other cases, a valid confidence interval may not exist.
mr_clr(mr_input(bx = calcium, bxse = calciumse,
by = fastgluc, byse = fastglucse, correl = calc.rho), nx=6351, ny=133010)##
## --------------------------------------------------------------------------
##
## 95 % confidence intervals using identification-robust methods
##
## Anderson--Rubin: [ 0.16 , 5.30 ]
## Kleibergen: [ 0.98 , 3.88 ]
## Conditional likelihood ratio: [ 0.96 , 3.92 ]
## --------------------------------------------------------------------------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.
MRAllObject_all <- mr_allmethods(MRInputObject, method = "all")
MRAllObject_all## Method Estimate Std Error 95% CI P-value
## Simple median 1.755 0.740 0.305 3.205 0.018
## Weighted median 2.683 0.419 1.862 3.504 0.000
## Penalized weighted median 2.681 0.420 1.857 3.505 0.000
##
## IVW 2.834 0.530 1.796 3.873 0.000
## [ reached 'max' / getOption("max.print") -- omitted 12 rows ]
MRAllObject_egger <- mr_allmethods(MRInputObject, method = "egger")
MRAllObject_egger## Method Estimate Std Error 95% CI P-value
## MR-Egger 3.253 0.770 1.743 4.762 0.000
## (intercept) -0.011 0.015 -0.041 0.018 0.451
## Penalized MR-Egger 3.421 0.531 2.380 4.461 0.000
## (intercept) -0.022 0.011 -0.043 0.000 0.051
## Robust MR-Egger 3.256 0.624 2.033 4.479 0.000
## [ reached 'max' / getOption("max.print") -- omitted 3 rows ]
MRAllObject_main <- mr_allmethods(MRInputObject, method = "main")
MRAllObject_main## Method Estimate Std Error 95% CI P-value
## Simple median 1.755 0.740 0.305 3.205 0.018
## Weighted median 2.683 0.419 1.862 3.504 0.000
## IVW 2.834 0.530 1.796 3.873 0.000
## MR-Egger 3.253 0.770 1.743 4.762 0.000
## (intercept) -0.011 0.015 -0.041 0.018 0.451Final 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.