library(MendelianRandomization)
## Warning: package 'MendelianRandomization' was built under R version 3.4.3
path.noproxy <- system.file("extdata", "vitD_snps_PhenoScanner.csv",
package = "MendelianRandomization")
path.proxies <- system.file("extdata", "vitD_snps_PhenoScanner_proxies.csv",
package = "MendelianRandomization")
# these two files from PhenoScanner are provided
# as part of the MendelianRandomization package
extract.pheno.csv(
exposure = "log(eGFR creatinine)", pmidE = 26831199, ancestryE = "European",
outcome = "Tanner stage", pmidO = 24770850, ancestryO = "European",
file = path.noproxy)
## SNP log(eGFR creatinine).beta log(eGFR creatinine).se
## 1 rs10741657 -0.0018 0.00092
## 2 rs17217119 -0.0066 0.00110
## 3 rs2282679 -0.0016 0.00100
## Tanner stage.beta Tanner stage.se
## 1 -0.0136 0.0135
## 2 0.0227 0.0165
## 3 -0.0042 0.0148
extract.pheno.csv(
exposure = "log(eGFR creatinine)", pmidE = 26831199, ancestryE = "European",
outcome = "Tanner stage", pmidO = 24770850, ancestryO = "European",
rsq.proxy = 0.6, file = path.proxies)
## Variant rs4944958 used as a proxy for rs12785878 in association with the outcome. R^2 value: 1.000.
## SNP log(eGFR creatinine).beta log(eGFR creatinine).se
## 1 rs10741657 -0.0018 0.00092
## 2 rs12785878 0.0024 0.00100
## 3 rs17217119 -0.0066 0.00110
## 4 rs2282679 -0.0016 0.00100
## Tanner stage.beta Tanner stage.se
## 1 -0.0136 0.0135
## 2 0.0102 0.0155
## 3 0.0227 0.0165
## 4 -0.0042 0.0148
extract.pheno.csv(
exposure = "log(eGFR creatinine)", pmidE = 26831199, ancestryE = "European",
outcome = "Asthma", pmidO = 20860503, ancestryO = "European",
rsq.proxy = 0.6, file = path.proxies)
## Variant rs2060793 used as a proxy for rs10741657 in association with the outcome. R^2 value: 1.000.
## Variant rs7944926 used as a proxy for rs12785878 in association with the outcome. R^2 value: 1.000.
## SNP log(eGFR creatinine).beta log(eGFR creatinine).se Asthma.beta
## 1 rs10741657 -0.0018 0.00092 0.00897
## 2 rs12785878 0.0024 0.00100 -0.00353
## 3 rs17217119 -0.0066 0.00110 -0.03906
## 4 rs2282679 -0.0016 0.00100 -0.01241
## Asthma.se
## 1 0.0202
## 2 0.0229
## 3 0.0255
## 4 0.0220
mr_allmethods(mr_input(bx = ldlc, bxse = ldlcse,
by = chdlodds, byse = chdloddsse), method="main", iterations = 100)
# iterations is set to 100 to reduce runtime for the mr_median method,
# at least 10000 iterations are recommended in practice
#MR-Egger method
mr_egger(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse))
##
## 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%
mr_egger(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse),
robust = TRUE)
##
## MR-Egger method
## (variants uncorrelated, random-effect model)
##
## Number of Variants = 28
## Robust model used.
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value
## MR-Egger 3.256 0.624 2.033, 4.479 0.000
## (intercept) -0.015 0.021 -0.055, 0.026 0.474
## ------------------------------------------------------------------
## Residual Standard Error : 1.498
## Heterogeneity test statistic = 97.7380 on 26 degrees of freedom, (p-value = 0.0000)
## I^2_GX statistic: 91.9%
mr_egger(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse),
penalized = TRUE)
##
## MR-Egger method
## (variants uncorrelated, random-effect model)
##
## Number of Variants = 28
## Weights of genetic variants with heterogeneous causal estimates have been penalized.
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value
## MR-Egger 3.421 0.531 2.380, 4.461 0.000
## (intercept) -0.022 0.011 -0.043, 0.000 0.051
## ------------------------------------------------------------------
## Residual Standard Error : 1.266
## Heterogeneity not calculated when weights are penalized.
## I^2_GX statistic: 91.9%
mr_egger(mr_input(calcium, calciumse, fastgluc, fastglucse, corr=calc.rho))
##
## MR-Egger method
## (variants correlated, random-effect model)
##
## Number of Variants = 6
##
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value
## MR-Egger 1.302 2.413 -3.427, 6.031 0.589
## (intercept) 0.009 0.021 -0.032, 0.050 0.666
## ------------------------------------------------------------------
## 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)
## correlated variants
#Inverse-variance weighted method
mr_ivw(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse))
##
## 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 = 99.5304 on 27 degrees of freedom, (p-value = 0.0000)
mr_ivw(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse),
robust = TRUE)
##
## Inverse-variance weighted method
## (variants uncorrelated, random-effect model)
##
## Number of Variants : 28
## Robust regression used.
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value
## IVW 2.797 0.307 2.195, 3.399 0.000
## ------------------------------------------------------------------
## Residual standard error = 1.987
## Heterogeneity test statistic = 106.5638 on 27 degrees of freedom, (p-value = 0.0000)
mr_ivw(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse),
penalized = TRUE)
##
## Inverse-variance weighted method
## (variants uncorrelated, random-effect model)
##
## Number of Variants : 28
## Weights of genetic variants with heterogeneous causal estimates have been penalized.
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value
## IVW 2.561 0.413 1.752, 3.370 0.000
## ------------------------------------------------------------------
## Residual standard error = 1.400
## Heterogeneity is not calculated when weights are penalized, or when there is only one variant in the analysis.
mr_ivw(mr_input(calcium, calciumse, fastgluc, fastglucse, corr=calc.rho))
## Warning in thetaIVW * Bx: Recycling array of length 1 in array-vector arithmetic is deprecated.
## Use c() or as.vector() instead.
##
## 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 = 2.0530 on 5 degrees of freedom, (p-value = 0.8418)
## correlated variants
#Maximum-likelihood method
mr_maxlik(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse))
##
## 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)
mr_maxlik(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse), psi=0.2)
##
## Maximum-likelihood method
## (variants uncorrelated, random-effect model)
##
## Number of Variants : 28
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value
## MaxLik 3.003 0.570 1.887, 4.120 0.000
## ------------------------------------------------------------------
## Residual standard error = 1.931
## Heterogeneity test statistic = 100.7226 on 27 degrees of freedom, (p-value = 0.0000)
mr_maxlik(mr_input(calcium, calciumse, fastgluc, fastglucse, corr=calc.rho))
##
## 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)
## correlated variants
#Median-based method
mr_median(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse),
weighting = "weighted", iterations = 100)
##
## Weighted median method
##
## Number of Variants : 28
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value
## Weighted median method 2.683 0.429 1.841, 3.524 0.000
## ------------------------------------------------------------------
# iterations is set to 100 to reduce runtime for the mr_median method,
# 10000 iterations are recommended in practice
mr_median(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse),
weighting = "simple", iterations = 100)
##
## Simple median method
##
## Number of Variants : 28
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value
## Simple median method 1.755 0.751 0.282, 3.228 0.020
## ------------------------------------------------------------------
mr_median(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse),
weighting = "penalized", iterations = 100)
##
## Penalized median method
##
## Number of Variants : 28
## ------------------------------------------------------------------
## Method Estimate Std Error 95% CI p-value
## Penalized median method 2.681 0.431 1.837, 3.525 0.000
## ------------------------------------------------------------------
#scatter plots of the genetic associations and/or causal estimates
mr_plot(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse),
line="egger", orientate = TRUE)
## We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
mr_plot(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse),
line="ivw", interactive=FALSE) # produces a static graph
mr_plot(mr_allmethods(mr_input(bx = ldlc, bxse = ldlcse,
by = chdlodds, byse = chdloddsse), method="all", iterations = 100))
# iterations is set to 100 to reduce runtime for the mr_median method,
# 10000 iterations are recommended in practice