TwoSampleMR 0.6.13

.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: >

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", , "[email protected]", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0003-0920-1055")),

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)

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)

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
+  ))
+}

R/mr.R

@@ -118,39 +118,39 @@ mr_method_list <- function()
 		list(
 			obj="mr_simple_median",
 			name="Simple median",
-			PubmedID="",
+			PubmedID="27061298",
 			Description="",
 			use_by_default=FALSE,
 			heterogeneity_test=FALSE
 		),
 		list(
 			obj="mr_weighted_median",
 			name="Weighted median",
-			PubmedID="",
+			PubmedID="27061298",
 			Description="",
 			use_by_default=TRUE,
 			heterogeneity_test=FALSE
 		),
 		list(
 			obj="mr_penalised_weighted_median",
 			name="Penalised weighted median",
-			PubmedID="",
+			PubmedID="27061298",
 			Description="",
 			use_by_default=FALSE,
 			heterogeneity_test=FALSE
 		),
 		list(
 			obj="mr_ivw",
 			name="Inverse variance weighted",
-			PubmedID="",
+			PubmedID="24114802",
 			Description="",
 			use_by_default=TRUE,
 			heterogeneity_test=TRUE
 		),
 		list(
 			obj = "mr_ivw_radial",
 			name = "IVW radial",
-			PubmedID = "",
+			PubmedID = "29961852",
 			Description = "",
 			use_by_default = FALSE,
 			heterogeneity_test = TRUE
@@ -174,31 +174,31 @@ mr_method_list <- function()
 		list(
 			obj="mr_simple_mode",
 			name="Simple mode",
-			PubmedID="",
+			PubmedID="29040600",
 			Description="",
 			use_by_default=TRUE,
 			heterogeneity_test=FALSE
 		),
 		list(
 			obj="mr_weighted_mode",
 			name="Weighted mode",
-			PubmedID="",
+			PubmedID="29040600",
 			Description="",
 			use_by_default=TRUE,
 			heterogeneity_test=FALSE
 		),
 		list(
 			obj="mr_weighted_mode_nome",
 			name="Weighted mode (NOME)",
-			PubmedID="",
+			PubmedID="29040600",
 			Description="",
 			use_by_default=FALSE,
 			heterogeneity_test=FALSE
 		),
 		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)

R/query.R

@@ -279,7 +279,8 @@
 	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))
 	{

R/read_data.R

@@ -165,12 +165,11 @@
 {
 
 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 @@
 		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))

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)))
+})

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)
+})

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: