Skip to contents

Mendelian Randomization wrapper (IVW, Egger, Weighted Median, Penalised WM)

Source: R/mr.R
mr.Rd

Mendelian Randomization wrapper (IVW, Egger, Weighted Median, Penalised WM)

Usage

mr(data, X, Y, alpha = 0.05, other_plots = FALSE)

Arguments

data

Data frame containing SNP summary statistics in wide format.

X

Character string. Exposure trait name.

Y

Character vector. Outcome trait names.

alpha

Significance level used for confidence intervals (default 0.05).

other_plots

Logical. If TRUE, produces metafor funnel and forest plots.

Value

Invisible list with components:

  • r Matrix of MR estimates for each exposure–outcome pair.

  • plots List of ggplot objects (MR scatter plots).

Details

Performs two–sample Mendelian Randomization using summary statistics stored in a wide data frame with columns named:

  • SNP

  • b.trait : effect estimate

  • SE.trait : standard error

For each outcome in Y, the function:

  1. Extracts instruments for exposure X

  2. Computes IVW, MR-Egger, Weighted Median and Penalised Weighted Median estimates

  3. Computes Cochran’s Q heterogeneity statistic

  4. Produces MR scatter plots

Methods implemented:

  • IVW (inverse-variance weighted regression)

  • MR-Egger regression

  • Weighted median estimator

  • Penalised weighted median

  • Cochran’s Q heterogeneity statistic (metafor)

Required packages: metafor, ggplot2, cowplot. Functions weighted.median() and mr.boot() must be available in the environment (e.g. from Mendelian randomization toolkits).

Examples

if (FALSE) { # \dontrun{
if (requireNamespace("metafor", quietly = TRUE) &&
    requireNamespace("ggplot2", quietly = TRUE) &&
    requireNamespace("cowplot", quietly = TRUE)) {
txt <- "
rs188743906  0.6804 0.1104  0.00177 0.01660        NA        NA
rs2289779   -0.0788 0.0134  0.00104 0.00261 -0.007543 0.0092258
rs117804300 -0.2281 0.0390 -0.00392 0.00855  0.109372 0.0362219
rs7033492   -0.0968 0.0147 -0.00585 0.00269  0.022793 0.0119903
rs10793962   0.2098 0.0212  0.00378 0.00536 -0.014567 0.0138196
rs635634    -0.2885 0.0153 -0.02040 0.00334  0.077157 0.0117123
rs176690    -0.0973 0.0142  0.00293 0.00306 -0.000007 0.0107781
rs147278971 -0.2336 0.0378 -0.01240 0.00792  0.079873 0.0397491
rs11562629   0.1155 0.0181  0.00960 0.00378 -0.010040 0.0151460
"
v <- c("SNP","b.LIF.R","SE.LIF.R","b.FEV1","SE.FEV1","b.CAD","SE.CAD")
mrdat <- setNames(as.data.frame(scan(text = txt,
                  what = list("",0,0,0,0,0,0), quiet=TRUE)), v)
res <- mr(mrdat, "LIF.R", c("CAD","FEV1"), other_plots=TRUE)
r <- as.data.frame(res$r, stringsAsFactors=FALSE)
rownames(r) <- r$IV
r$IV <- NULL
r[] <- lapply(r, function(x) as.numeric(as.character(x)))
b_cols  <- grep("^b[A-Z]", names(r), value=TRUE)
se_cols <- paste0("se", substring(b_cols, 2))
keep <- se_cols %in% names(r)
b_cols <- b_cols[keep]; se_cols <- se_cols[keep]
if(length(b_cols) > 0){
  pvals <- sapply(seq_along(b_cols), function(i)
    2 * pnorm(-abs(r[[b_cols[i]]] / r[[se_cols[i]]])))
  pvals <- as.data.frame(pvals)
  colnames(pvals) <- substring(b_cols, 2)
  pvals[] <- lapply(pvals, format.pval, digits=3, eps=1e-4)
  results_table <- cbind(round(r,3), pvals)
} else {
  results_table <- round(r,3)
}
results_table
res$plots
}
} # }