diff --git a/.github/workflows/check-full.yaml b/.github/workflows/check-full.yaml index c8dea762..64317ec5 100644 --- a/.github/workflows/check-full.yaml +++ b/.github/workflows/check-full.yaml @@ -31,9 +31,9 @@ jobs: - {os: macos-latest, r: 'release'} - {os: macos-13, r: 'release'} - {os: windows-latest, r: 'release'} - - {os: ubuntu-latest, r: 'devel', http-user-agent: 'release'} - {os: ubuntu-latest, r: 'release'} - {os: ubuntu-latest, r: 'oldrel-1'} + - {os: ubuntu-latest, r: 'devel', http-user-agent: 'release'} - {os: ubuntu-latest, r: '4.3.2'} - {os: ubuntu-24.04-arm, r: 'release', rspm: 'no' } @@ -53,6 +53,10 @@ jobs: http-user-agent: ${{ matrix.config.http-user-agent }} use-public-rspm: ${{ matrix.config.rspm || 'true' }} + - name: Set pak options + shell: bash + run: echo 'options(pkg.sysreqs_db_update_timeout = as.difftime(59, units = "secs"))' >> ~/.Rprofile + - uses: r-lib/actions/setup-r-dependencies@v2 with: extra-packages: > diff --git a/DESCRIPTION b/DESCRIPTION index 014bcde4..2b7c824a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: TwoSampleMR Title: Two Sample MR Functions and Interface to MRC Integrative Epidemiology Unit OpenGWAS Database -Version: 0.6.12 +Version: 0.6.13 Authors@R: c( person("Gibran", "Hemani", , "g.hemani@bristol.ac.uk", role = c("aut", "cre"), comment = c(ORCID = "0000-0003-0920-1055")), diff --git a/NAMESPACE b/NAMESPACE index 7e7dedac..4e35877c 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -51,6 +51,7 @@ export(mr_egger_regression) export(mr_egger_regression_bootstrap) export(mr_forest_plot) export(mr_funnel_plot) +export(mr_grip) export(mr_heterogeneity) export(mr_ivw) export(mr_ivw_fe) diff --git a/NEWS.md b/NEWS.md index b32ec49b..f3ad7d97 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,12 @@ +# TwoSampleMR v0.6.13 + +(Release date 2025-03-26) + +* Added `mr_grip()` function which implements the MR-GRIP (modified MR-Egger with the Genotype Recoding Invariance Property) method of Dudbridge and Bowden et al. (2025). + The new method can be accessed by `mr(dat, method_list = "mr_grip")` or it can be added to the default list of methods with `mr(dat, method_list = c(subset(mr_method_list(), use_by_default)$obj, "mr_grip"))`. +* Added Pub Med IDs for more of the methods. +* The `format_data()` function no longer causes a stack overflow when its `dat` argument is not a variable (thanks to @DarwinAwardWinner) + # TwoSampleMR v0.6.12 (Release date 2025-03-18) diff --git a/R/mr-grip.R b/R/mr-grip.R new file mode 100644 index 00000000..3b5f9cbb --- /dev/null +++ b/R/mr-grip.R @@ -0,0 +1,79 @@ +#' MR-GRIP: a modified MR-Egger model with the Genotype Recoding Invariant Property +#' +#' This implements the modified MR-Egger model with the Genotype Recoding Invariant Property (MR-GRIP) due to Dudbridge and Bowden et al. (2025). +#' It is well known that the results of MR-Egger are sensitive to which alleles are designated as the effect alleles. +#' A pragmatic convention is to orient all SNPs to have positive effects on the exposure, which has some advantages in interpretation but also brings some philosophical limitations. +#' The MR-GRIP model is a modification to the MR-Egger model in which each term is multiplied by the genotype-phenotype associations. +#' This makes each term in the model invariant to allele coding. +#' +#' @param b_exp Vector of genetic effects on exposure. +#' @param b_out Vector of genetic effects on outcome. +#' @param se_exp Standard errors of genetic effects on exposure. +#' @param se_out Standard errors of genetic effects on outcome. +#' @param parameters List of parameters. +#' +#' @export +#' @return List of with the following elements: +#' \describe{ +#' \item{b}{MR estimate} +#' \item{se}{Standard error of MR estimate} +#' \item{pval}{p-value of MR estimate} +#' \item{Q, Q_df, Q_pval}{Heterogeneity stats} +#' \item{b.wi}{MR estimate adjusting for weak instruments} +#' \item{se.wi}{Standard error adjusting for weak instruments} +#' \item{pval.wi}{p-value adjusting for weak instruments} +#' \item{mod}{Summary of regression} +#' \item{dat}{Original data used for MR-GRIP} +#' } +mr_grip <- function(b_exp, b_out, se_exp, se_out, parameters) { + if (length(b_exp) != length(b_out)) stop("The lengths of b_exp and b_out are not equal.") + if (length(se_exp) != length(se_out)) stop("The lengths of se_exp and se_out are not equal.") + if (length(b_exp) != length(se_out)) stop("The lengths of b_exp and se_out are not equal.") + + nulllist <- list( + b = NA, + se = NA, + pval = NA, + nsnp = NA, + Q = NA, + Q_df = NA, + Q_pval = NA, + mod = NA, + smod = NA, + dat = NA + ) + if ( + sum(!is.na(b_exp) & !is.na(b_out) & !is.na(se_exp) & !is.na(se_out)) < 3 + ) { + return(nulllist) + } + + dat <- data.frame( + b_out = b_out, + b_exp = b_exp, + se_exp = se_exp, + se_out = se_out + ) + grip_out <- b_out * b_exp + grip_exp <- b_exp^2 + # GRIP regression. Includes intercept. Weights designed to replicate IVW under no intercept. + mod <- stats::lm(grip_out ~ grip_exp, weights = 1 / (grip_exp * se_out^2)) + smod <- summary(mod) + b <- stats::coefficients(smod)[2, 1] + se <- stats::coefficients(smod)[2, 2] + b.adj <- NA + se.adj <- NA + pval.adj <- NA + pval <- 2 * stats::pt(abs(b / se), length(b_exp) - 2L, lower.tail = FALSE) + return(list( + b = b, + se = se, + pval = pval, + b.adj = b.adj, + se.adj = se.adj, + pval.adj = pval.adj, + nsnp = length(b_exp), + mod = smod, + dat = dat + )) +} diff --git a/R/mr.R b/R/mr.R index 657d1432..3691142d 100644 --- a/R/mr.R +++ b/R/mr.R @@ -118,7 +118,7 @@ mr_method_list <- function() list( obj="mr_simple_median", name="Simple median", - PubmedID="", + PubmedID="27061298", Description="", use_by_default=FALSE, heterogeneity_test=FALSE @@ -126,7 +126,7 @@ mr_method_list <- function() list( obj="mr_weighted_median", name="Weighted median", - PubmedID="", + PubmedID="27061298", Description="", use_by_default=TRUE, heterogeneity_test=FALSE @@ -134,7 +134,7 @@ mr_method_list <- function() list( obj="mr_penalised_weighted_median", name="Penalised weighted median", - PubmedID="", + PubmedID="27061298", Description="", use_by_default=FALSE, heterogeneity_test=FALSE @@ -142,7 +142,7 @@ mr_method_list <- function() list( obj="mr_ivw", name="Inverse variance weighted", - PubmedID="", + PubmedID="24114802", Description="", use_by_default=TRUE, heterogeneity_test=TRUE @@ -150,7 +150,7 @@ mr_method_list <- function() list( obj = "mr_ivw_radial", name = "IVW radial", - PubmedID = "", + PubmedID = "29961852", Description = "", use_by_default = FALSE, heterogeneity_test = TRUE @@ -174,7 +174,7 @@ mr_method_list <- function() list( obj="mr_simple_mode", name="Simple mode", - PubmedID="", + PubmedID="29040600", Description="", use_by_default=TRUE, heterogeneity_test=FALSE @@ -182,7 +182,7 @@ mr_method_list <- function() list( obj="mr_weighted_mode", name="Weighted mode", - PubmedID="", + PubmedID="29040600", Description="", use_by_default=TRUE, heterogeneity_test=FALSE @@ -190,7 +190,7 @@ mr_method_list <- function() list( obj="mr_weighted_mode_nome", name="Weighted mode (NOME)", - PubmedID="", + PubmedID="29040600", Description="", use_by_default=FALSE, heterogeneity_test=FALSE @@ -198,7 +198,7 @@ mr_method_list <- function() list( obj="mr_simple_mode_nome", name="Simple mode (NOME)", - PubmedID="", + PubmedID="29040600", Description="", use_by_default=FALSE, heterogeneity_test=FALSE @@ -226,6 +226,14 @@ mr_method_list <- function() Description="Doesn't use any weights", use_by_default=FALSE, heterogeneity_test=TRUE + ), + list( + obj = "mr_grip", + name = "MR GRIP", + PubmedID = "", + Description = "Allele coding invariant regression", + use_by_default = FALSE, + heterogeneity_test = FALSE ) ) a <- lapply(a, as.data.frame) diff --git a/R/query.R b/R/query.R index 45006818..be8f80e2 100644 --- a/R/query.R +++ b/R/query.R @@ -279,7 +279,8 @@ format_d <- function(d) d$mr_keep.outcome <- apply(d[, mrcols], 1, function(x) !any(is.na(x))) if(any(!d$mr_keep.outcome)) { - warning("The following SNP(s) are missing required information for the MR tests and will be excluded\n", paste(subset(d, !mr_keep.outcome)$SNP, collapse="\n")) + missinginfosnps <- paste(subset(d, !mr_keep.outcome)$SNP, collapse = " ") + warning("The following SNP(s) are missing required information for the MR tests and will be excluded: ", missinginfosnps) } if(all(!d$mr_keep.outcome)) { diff --git a/R/read_data.R b/R/read_data.R index 8741f859..b0bbc250 100644 --- a/R/read_data.R +++ b/R/read_data.R @@ -165,12 +165,11 @@ format_data <- function(dat, type="exposure", snps=NULL, header=TRUE, { if (inherits(dat, "data.table")) { - datname <- deparse(substitute(dat)) - stop(paste0( - "Your ", datname, " data.frame is also of class 'data.table', ", - "please reformat as simply a data.frame with ", datname, " <- data.frame(", - datname, ") and then rerun your format_data() call." - )) + stop( + "Your data.frame is also of class 'data.table' ", + "please reformat as simply a data.frame with data.frame() or as.data.frame() ", + "and then rerun format_data()." + ) } all_cols <- c(phenotype_col, snp_col, beta_col, se_col, eaf_col, effect_allele_col, other_allele_col, pval_col, units_col, ncase_col, ncontrol_col, samplesize_col, gene_col, id_col, z_col, info_col, chr_col, pos_col) @@ -482,7 +481,8 @@ if (inherits(dat, "data.table")) { dat$mr_keep.outcome <- dat$mr_keep.outcome & apply(dat[, mrcols_present], 1, function(x) !any(is.na(x))) if(any(!dat$mr_keep.outcome)) { - warning("The following SNP(s) are missing required information for the MR tests and will be excluded\n", paste(subset(dat, !mr_keep.outcome)$SNP, collapse="\n")) + missinginfosnps <- paste(subset(dat, !mr_keep.outcome)$SNP, collapse = " ") + warning("The following SNP(s) are missing required information for the MR tests and will be excluded: ", missinginfosnps) } } if(all(!dat$mr_keep.outcome)) diff --git a/man/mr_grip.Rd b/man/mr_grip.Rd new file mode 100644 index 00000000..b5b3f235 --- /dev/null +++ b/man/mr_grip.Rd @@ -0,0 +1,40 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/mr-grip.R +\name{mr_grip} +\alias{mr_grip} +\title{MR-GRIP: a modified MR-Egger model with the Genotype Recoding Invariant Property} +\usage{ +mr_grip(b_exp, b_out, se_exp, se_out, parameters) +} +\arguments{ +\item{b_exp}{Vector of genetic effects on exposure.} + +\item{b_out}{Vector of genetic effects on outcome.} + +\item{se_exp}{Standard errors of genetic effects on exposure.} + +\item{se_out}{Standard errors of genetic effects on outcome.} + +\item{parameters}{List of parameters.} +} +\value{ +List of with the following elements: +\describe{ +\item{b}{MR estimate} +\item{se}{Standard error of MR estimate} +\item{pval}{p-value of MR estimate} +\item{Q, Q_df, Q_pval}{Heterogeneity stats} +\item{b.wi}{MR estimate adjusting for weak instruments} +\item{se.wi}{Standard error adjusting for weak instruments} +\item{pval.wi}{p-value adjusting for weak instruments} +\item{mod}{Summary of regression} +\item{dat}{Original data used for MR-GRIP} +} +} +\description{ +This implements the modified MR-Egger model with the Genotype Recoding Invariant Property (MR-GRIP) due to Dudbridge and Bowden et al. (2025). +It is well known that the results of MR-Egger are sensitive to which alleles are designated as the effect alleles. +A pragmatic convention is to orient all SNPs to have positive effects on the exposure, which has some advantages in interpretation but also brings some philosophical limitations. +The MR-GRIP model is a modification to the MR-Egger model in which each term is multiplied by the genotype-phenotype associations. +This makes each term in the model invariant to allele coding. +} diff --git a/tests/testthat/test_format_data.R b/tests/testthat/test_format_data.R index d67b85f4..b2ae2b4e 100644 --- a/tests/testthat/test_format_data.R +++ b/tests/testthat/test_format_data.R @@ -45,3 +45,8 @@ test_that("format_data() should not error after having its data.table class remo samplesize_col = "n" )) }) + +test_that("format_data() should not cause a stack overflow", { + a <- data.table::data.table(x = sample(1:1e6)) + expect_error(do.call(format_data, list(a))) +}) diff --git a/tests/testthat/test_mr.R b/tests/testthat/test_mr.R index c07d259d..ce33ec96 100644 --- a/tests/testthat/test_mr.R +++ b/tests/testthat/test_mr.R @@ -51,3 +51,10 @@ test_that("mr.raps shrinkage option", { expect_equal(ncol(res5), 9L) expect_equal(res5[1, "b"], 0.4647, tolerance = 1e-3) }) + +test_that("mr_grip()", { + res6 <- suppressWarnings(mr(dat, method_list = "mr_grip")) + expect_equal(nrow(res6), 1L) + expect_equal(ncol(res6), 9L) + expect_equal(res6[1, "b"], 0.490, tolerance = 1e-3) +}) diff --git a/vignettes/introduction.Rmd b/vignettes/introduction.Rmd index 84a6dfb3..2ca16a4d 100644 --- a/vignettes/introduction.Rmd +++ b/vignettes/introduction.Rmd @@ -32,8 +32,6 @@ The general principles [@DaveySmith2003; @DaveySmithHemani2014], and statistical This package uses the [ieugwasr](https://github.com/mrcieu/ieugwasr) package to connect to the database of thousands of complete GWAS summary data. -* * * - ## Installation To install directly from the GitHub repository do the following: @@ -45,8 +43,6 @@ install_github("MRCIEU/TwoSampleMR") If you don't have the `remotes` package install it from CRAN using `install.packages("remotes")`. -* * * - ## Overview The workflow for performing MR is as follows: diff --git a/vignettes/perform_mr.Rmd b/vignettes/perform_mr.Rmd index 0a17c48d..080f136b 100644 --- a/vignettes/perform_mr.Rmd +++ b/vignettes/perform_mr.Rmd @@ -84,8 +84,6 @@ mr(dat, method_list = c("mr_egger_regression", "mr_ivw")) By default, all the methods that are labelled `TRUE` in the `use_by_default` column are used by the `mr()` function. -* * * - ## Sensitivity analyses ### Heterogeneity statistics @@ -144,8 +142,6 @@ res_loo <- mr_leaveoneout(dat) By default the method used is the inverse variance weighted method, but this can be changed by using the `method` argument. -* * * - ## Plots There are a few ways to visualise the results, listed below @@ -478,18 +474,16 @@ plot1 dev.off() ``` -* * * +## MR-RAPS: Many weak instruments analysis -## MR.RAPS: Many weak instruments analysis +MR-RAPS (Robust Adjusted Profile Score) is a recently proposed method that considers the measurement error in SNP-exposure effects, is unbiased when there are many (e.g. hundreds of) weak instruments, and is robust to systematic and idiosyncratic pleiotropy. See @zhao2020 for more detail about the statistical methodology. -MR.RAPS (Robust Adjusted Profile Score) is a recently proposed method that considers the measurement error in SNP-exposure effects, is unbiased when there are many (e.g. hundreds of) weak instruments, and is robust to systematic and idiosyncratic pleiotropy. See @zhao2020 for more detail about the statistical methodology. - -MR.RAPS is implemented in the R package **mr.raps** that is available on [GitHub](https://github.com/qingyuanzhao/mr.raps) which is installed when you install TwoSampleMR. It can be directly called from TwoSampleMR by +MR-RAPS is implemented in the R package **mr.raps** that is available on [GitHub](https://github.com/qingyuanzhao/mr.raps) which is installed when you install TwoSampleMR. It can be directly called from TwoSampleMR by ```{r eval=FALSE} res <- mr(dat, method_list = c("mr_raps")) ``` -MR.RAPS comes with two main options: `over.dispersion` (whether the method should consider systematic pleiotropy) and `loss.function` (either `"l2"`, `"huber"`, or `"tukey"`). The latter two loss functions are robust to idiosyncratic pleiotropy. The default option is `over.dispersion = TRUE` and `loss.function = "huber"`. To change these options, modify the `parameters` argument of `mr()` by (for example) +MR-RAPS comes with two main options: `over.dispersion` (whether the method should consider systematic pleiotropy) and `loss.function` (either `"l2"`, `"huber"`, or `"tukey"`). The latter two loss functions are robust to idiosyncratic pleiotropy. The default option is `over.dispersion = TRUE` and `loss.function = "huber"`. To change these options, modify the `parameters` argument of `mr()` by (for example) ```{r eval=FALSE} res <- mr( @@ -499,7 +493,19 @@ res <- ) ``` -* * * +## MR-GRIP + +MR-GRIP is a recently proposed method due to Dudbridge and Bowden et al. (2025) which is a modification to the MR-Egger method which has the Genotype Recoding Invariance Property (GRIP), it can be implemented using the code below. + +```{r eval=FALSE} +res <- mr(dat, method_list = c("mr_grip")) +``` + +Or include it with the default set of methods as follows. + +```{r, eval=FALSE} +mr(dat, method_list = c(subset(mr_method_list(), use_by_default)$obj, "mr_grip")) +``` ## Reports @@ -513,8 +519,6 @@ By default this produces a html file in the current working directory, but see t This function will create a separate report file for every exposure-outcome combination that is present in the `dat` object. -* * * - ## MR Steiger directionality test This is an implementation of the method described here: @@ -549,8 +553,6 @@ mr_steiger( ) ``` -* * * - ## Multivariable MR When SNPs instrument multiple potential exposures, for example in the case of different lipid fractions, one method for overcoming this problem is to estimate the influence of each lipid conditioning on the effects of the SNPs on the other lipids. Multivariable MR can be performed using the R package as follows. Practically speaking, this is the process that needs to occur from the perspective of generating the data in the correct format: @@ -615,8 +617,6 @@ The current plots being generated are not necessarily adequate because while the If you want to perform analysis with your local summary data (i.e. not in the OpenGWAS database) then use then look up the `mv_extract_exposures_local()` function instead of the `mv_extract_exposures()` function. -* * * - ## MR estimates when instruments are correlated In the examples shown so far it is assumed that instruments are independent (i.e. they are not in linkage disequilibrium, LD). This is to avoid 'double counting' effects. An alternative approach is to estimate the MR effects accounting for the correlation between variants. @@ -660,8 +660,6 @@ dat2 <- try(dat_to_MRInput(dat, get_correlation = TRUE)) if (class(dat2) != "try-error") MendelianRandomization::mr_ivw(dat2[[1]], correl = TRUE) ``` -* * * - ## MR-MoE: Using a mixture of experts machine learning approach We recently developed MR-MoE, a method to choose the most appropriate amongst several MR tests using a machine learning algorithm. Note that the method is still under review, but full details are described here: . @@ -714,8 +712,6 @@ The result object itself is a list with the following elements: Looking at the `estimates`, we see that there is a column called `MOE` which is the predicted AUROC curve performance of each method. -* * * - ## Post MR results management The TwoSampleMR package also provides the following functions for managing or editing MR results.