TwoSampleMR 0.7.13

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: TwoSampleMR
 Title: Two Sample MR Functions and Interface to MRC Integrative
     Epidemiology Unit OpenGWAS Database
-Version: 0.7.12
+Version: 0.7.13
 Authors@R: c(
     person("Gibran", "Hemani", , "[email protected]", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0003-0920-1055")),
@@ -69,5 +69,5 @@ Remotes:
     qingyuanzhao/mr.raps,
     WSpiller/RadialMR
 Config/roxygen2/markdown: TRUE
-Config/roxygen2/version: 8.1.0
+Config/roxygen2/version: 8.1.1
 Encoding: UTF-8

NEWS.md

@@ -1,3 +1,16 @@
+# TwoSampleMR v0.7.13
+
+(Release date 2026-10-20)
+
+* Fixed URL in gwas2020 vignette
+* Clarified the steps in the Multivariable MR section of the Perform MR vignette and fixed the object passed to `extract_outcome_data()` in its example (thanks @hdraisma, #340)
+* Fixed `harmonise_ld_dat()`, and hence `dat_to_MRInput(get_correlations = TRUE)`, failing when the LD matrix SNPs were in a different order to the summary data; the LD matrix is now reordered to match (thanks @marwaalarab, #348)
+* Fixed small p-values being shown as 0 in the tables in `mr_report()` reports; p-values are now shown to 3 significant figures (thanks @sonia-shah, #42)
+* `run_mrmix()` now removes SNPs with an `eaf.exposure` that is missing or not strictly between 0 and 1, with a warning, instead of erroring in `MRMix::standardize()` (#169)
+* The mode based estimators (`mr_mode()`, `mr_simple_mode()`, `mr_weighted_mode()`, `mr_simple_mode_nome()`, and `mr_weighted_mode_nome()`) now remove SNPs with a `beta.exposure` of 0, with a warning, instead of erroring with "non-finite 'bw'" in `density()` (thanks @Julie1130, #570)
+* In `mv_multiple()` and `mv_residual()` with `instrument_specific = FALSE` (the default), the `nsnp` column now reports the number of SNPs used in the analysis, which is all the SNPs, rather than the number passing `pval_threshold` for each exposure. In `mv_residual()` this also fixes estimates being set to `NA` when few SNPs passed `pval_threshold` for an exposure, even though all SNPs were used (thanks @MaVdb and @mvab, #78)
+* Expanded the documentation of `effective_n()` to give its formula, explain that it is half the effective sample size used by METAL, and where it is used in the package (thanks @carbocation, #433)
+
 # TwoSampleMR v0.7.12
 
 (Release date 2026-10-08)
@@ -69,7 +82,7 @@
 
 (Release date 2026-05-04)
 
-* Amend <https://gwas.mrcieu.ac.uk> URLs to <https://opengwas.io/>
+* Amend `https://gwas.mrcieu.ac.uk` URLs to <https://opengwas.io/>
 * Bump version of roxygen2 to 8.0.0 and regenerate documentation
 * Replace `array(1:nexp)` placeholder allocations with `numeric(nexp)`/`integer(nexp)` in multivariable MR functions
 * Replace `1:n` with `seq_len(n)` in loop indices across `singlesnp`, `leaveoneout`, `rucker`, `mr_mode`, and `multivariable_mr` for empty-vector safety

R/add_rsq.r

@@ -379,7 +379,17 @@ get_population_allele_frequency <- function(af, prop, odds_ratio, prevalence) {
 
 #' Estimate the effective sample size in a case control study
 #'
-#' Taken from <https://www.nature.com/articles/nprot.2014.071>
+#' Computes the effective sample size as the harmonic mean of the number of cases and controls,
+#' \deqn{N_{eff} = \frac{2}{1/N_{cases} + 1/N_{controls}},}{N_eff = 2 / (1/N_cases + 1/N_controls),}
+#' as given by Winkler et al. (2014) <https://www.nature.com/articles/nprot.2014.071>.
+#'
+#' Note that other software, such as METAL (Willer et al., 2010, \doi{10.1093/bioinformatics/btq340}),
+#' uses \eqn{4/(1/N_{cases} + 1/N_{controls})}{4 / (1/N_cases + 1/N_controls)}, which is twice this value.
+#' With that definition a balanced study (equal numbers of cases and controls) has an effective sample size equal to its total sample size,
+#' whereas with this function it is equal to half the total sample size, i.e., the number of cases.
+#' If you require the METAL definition multiply the result of this function by 2.
+#'
+#' This function is used by [add_rsq()] to obtain the sample sizes for binary traits (those with units of `"log odds"`) which are then used by [directionality_test()] and [steiger_filtering()].
 #'
 #' @param ncase Vector of number of cases
 #' @param ncontrol Vector of number of controls

R/knit.R

@@ -150,6 +150,7 @@ mr_report <- function(
         }
       }
     )
+    tablist <- lapply(tablist, format_report_pvals)
 
     plotlist <- lapply(p, function(x) {
       d <- attributes(x)$split_labels
@@ -172,6 +173,19 @@ mr_report <- function(
   return(output_file)
 }
 
+# Format p-value columns as character so that kable() does not round small
+# p-values to 0
+format_report_pvals <- function(x) {
+  if (is.null(x)) {
+    return(NULL)
+  }
+  i <- grepl("pval", names(x)) & vapply(x, is.numeric, logical(1))
+  x[i] <- lapply(x[i], function(p) {
+    ifelse(is.na(p), NA_character_, trimws(formatC(p, format = "g", digits = 3)))
+  })
+  x
+}
+
 sanitise_string <- function(x) {
   gsub(" ", "_", gsub("[^[:alnum:] ]", "", x))
 }

R/mr_mode.R

@@ -13,9 +13,24 @@ mr_mode <- function(dat, parameters = default_parameters(), mode_method = "all")
     dat <- subset(dat, mr_keep)
   }
 
+  # The ratio estimate is infinite for SNPs with a beta.exposure of 0
+  index <- dat$beta.exposure == 0
+  if (any(index)) {
+    warning(
+      "Removing ",
+      sum(index),
+      " SNP(s) with beta.exposure of 0 from the mode based estimates"
+    )
+    dat <- dat[!index, ]
+  }
+
   if (nrow(dat) < 3) {
     warning("Need at least 3 SNPs")
-    return(NULL)
+    if (mode_method == "all") {
+      return(NULL)
+    } else {
+      return(list(b = NA, se = NA, pval = NA, nsnp = NA))
+    }
   }
 
   b_exp <- dat$beta.exposure

R/multivariable_mr.R

@@ -478,7 +478,7 @@ mv_harmonise_data <- function(exposure_dat, outcome_dat, harmonise_strictness =
 #' @param mvdat Output from [mv_harmonise_data()].
 #' @param intercept Should the intercept by estimated (`TRUE`) or force line through the origin (`FALSE`, default).
 #' @param instrument_specific Should the estimate for each exposure be obtained by using all instruments from all exposures (`FALSE`, default) or by using only the instruments specific to each exposure (`TRUE`).
-#' @param pval_threshold P-value threshold to include instruments. The default is `5e-8`.
+#' @param pval_threshold P-value threshold to include instruments. Only used when `instrument_specific = TRUE`. The default is `5e-8`.
 #' @param plots Create plots? The default is `FALSE`.
 #'
 #' @export
@@ -505,8 +505,13 @@ mv_residual <- function(
   nom <- colnames(beta.exposure)
   nom2 <- mvdat$expname$exposure[match(nom, mvdat$expname$id.exposure)]
   for (i in seq_len(nexp)) {
-    # For this exposure, only keep SNPs that meet some p-value threshold
-    index <- pval.exposure[, i] < pval_threshold
+    # For this exposure, only keep SNPs that meet some p-value threshold if the
+    # estimate is instrument specific, otherwise all SNPs are used
+    if (instrument_specific) {
+      index <- pval.exposure[, i] < pval_threshold
+    } else {
+      index <- rep(TRUE, nrow(beta.exposure))
+    }
 
     # Get outcome effects adjusted for all effects on all other exposures
     if (intercept) {
@@ -586,7 +591,7 @@ mv_residual <- function(
 #' @param mvdat Output from [mv_harmonise_data()].
 #' @param intercept Should the intercept by estimated (`TRUE`) or force line through the origin (`FALSE`, default).
 #' @param instrument_specific Should the estimate for each exposure be obtained by using all instruments from all exposures (`FALSE`, default) or by using only the instruments specific to each exposure (`TRUE`).
-#' @param pval_threshold P-value threshold to include instruments. The default is `5e-8`.
+#' @param pval_threshold P-value threshold to include instruments. Only used when `instrument_specific = TRUE`. The default is `5e-8`.
 #' @param plots Create plots? The default is `FALSE`.
 #'
 #' @export
@@ -614,8 +619,13 @@ mv_multiple <- function(
   nom <- colnames(beta.exposure)
   nom2 <- mvdat$expname$exposure[match(nom, mvdat$expname$id.exposure)]
   for (i in seq_len(nexp)) {
-    # For this exposure, only keep SNPs that meet some p-value threshold
-    index <- pval.exposure[, i] < pval_threshold
+    # For this exposure, only keep SNPs that meet some p-value threshold if the
+    # estimate is instrument specific, otherwise all SNPs are used
+    if (instrument_specific) {
+      index <- pval.exposure[, i] < pval_threshold
+    } else {
+      index <- rep(TRUE, nrow(beta.exposure))
+    }
 
     # # Get outcome effects adjusted for all effects on all other exposures
     # marginal_outcome[,i] <- lm(beta.outcome ~ beta.exposure[, -c(i)])$res

R/other_formats.R

@@ -78,10 +78,14 @@ dat_to_MRInput <- function(dat, get_correlations = FALSE, pop = "EUR") {
 harmonise_ld_dat <- function(x, ld) {
   snpnames <- do.call(rbind, strsplit(rownames(ld), split = "_"))
   i1 <- snpnames[, 1] %in% x$SNP
-  ld <- ld[i1, i1]
-  snpnames <- snpnames[i1, ]
+  ld <- ld[i1, i1, drop = FALSE]
+  snpnames <- snpnames[i1, , drop = FALSE]
   i2 <- x$SNP %in% snpnames[, 1]
   x <- x[i2, ]
+  # Put the LD matrix in the same SNP order as x
+  i3 <- match(x$SNP, snpnames[, 1])
+  ld <- ld[i3, i3, drop = FALSE]
+  snpnames <- snpnames[i3, , drop = FALSE]
   stopifnot(all(snpnames[, 1] == x$SNP))
   x$effect_allele.exposure <- as.character(x$effect_allele.exposure)
   x$other_allele.exposure <- as.character(x$other_allele.exposure)
@@ -125,7 +129,7 @@ harmonise_ld_dat <- function(x, ld) {
 
   if (!all(snpnames$keep)) {
     message("Removing ", sum(!snpnames$keep), " variants due to harmonisation issues")
-    ld <- ld[snpnames$keep, snpnames$keep]
+    ld <- ld[snpnames$keep, snpnames$keep, drop = FALSE]
     x <- x[snpnames$keep, ]
   }
   return(list(x = x, ld = ld))
@@ -251,6 +255,7 @@ mr_ivw_radial <- function(b_exp, b_out, se_exp, se_out, parameters = default_par
 #' See <https://github.com/gqi/MRMix> for more details.
 #'
 #' @param dat Output from [harmonise_data()]. Ensures that no eaf.exposure values are missing.
+#' Missing eaf.exposure values are replaced by eaf.outcome values, and SNPs with an eaf.exposure that is still missing or is not strictly between 0 and 1 are removed with a warning.
 #'
 #' @export
 #' @return List of results, with one list item for every exposure/outcome pair in dat object
@@ -283,6 +288,24 @@ run_mrmix <- function(dat) {
       ))
       x$eaf.exposure[index] <- x$eaf.outcome[index]
     }
+    # standardize() scales by sqrt(2 * eaf * (1 - eaf)), so SNPs with an eaf of
+    # 0 or 1 would give standard errors of 0
+    index <- is.na(x$eaf.exposure) | x$eaf.exposure <= 0 | x$eaf.exposure >= 1
+    if (any(index)) {
+      warning(paste0(
+        x$id.exposure[1],
+        ".",
+        x$id.outcome[1],
+        " - Removing ",
+        sum(index),
+        " SNP(s) with eaf.exposure missing or not strictly between 0 and 1"
+      ))
+      x <- x[!index, ]
+    }
+    if (nrow(x) == 0) {
+      warning(paste0(combos$id.exposure[i], ".", combos$id.outcome[i], " - No SNPs remaining"))
+      return(NULL)
+    }
     l <- MRMix::standardize(
       x$beta.exposure,
       x$beta.outcome,

inst/CITATION

@@ -30,7 +30,7 @@ bibentry(
   volume = "7",
   journal = "eLife",
   pages = "e34408",
-  url = "https://elifesciences.org/articles/34408",
+  url = "https://doi.org/10.7554/eLife.34408",
   doi = "10.7554/eLife.34408",
   key = "twosamplemr",
   textVersion = paste(

justfile

@@ -14,3 +14,5 @@ dev:
     Rscript -e "pak::local_install_dev_deps()"
 readme:
     Rscript -e "rmarkdown::render('README.Rmd', output_options = list(html_preview = FALSE))"
+urlcheck:
+    Rscript -e "urlchecker::url_check()"

tests/testthat/test_knit.R

@@ -0,0 +1,15 @@
+test_that("format_report_pvals keeps small p-values", {
+  x <- data.frame(
+    method = c("a", "b", "c"),
+    b = c(0.5, 0.1, 0.2),
+    pval = c(1e-50, 0.0123, NA),
+    Q_pval = c(0.5, 3.2e-12, 1)
+  )
+  out <- format_report_pvals(x)
+  expect_equal(out$pval, c("1e-50", "0.0123", NA))
+  expect_equal(out$Q_pval, c("0.5", "3.2e-12", "1"))
+  expect_equal(out$b, x$b)
+  expect_equal(out$method, x$method)
+  expect_null(format_report_pvals(NULL))
+  expect_false(grepl("\\| *0 *\\|", paste(knitr::kable(out), collapse = "\n")))
+})

tests/testthat/test_mr_mode.R

@@ -41,3 +41,37 @@ test_that("mr_mode() simple-mode SE is not inflated when nsnp divides nboot (iss
   # Correct simple-mode SE is ~0.17; the bug inflated it to ~0.57.
   expect_lt(m$se[m$method == "Simple mode"], 0.35)
 })
+
+test_that("mr_mode() removes SNPs with beta.exposure of 0 (issue #570)", {
+  d <- dat
+  d$beta.exposure[1:2] <- 0
+  set.seed(1234)
+  expect_warning(m <- mr_mode(d), "Removing 2 SNP\\(s\\) with beta.exposure of 0")
+  set.seed(1234)
+  m0 <- mr_mode(dat[-(1:2), ])
+  expect_equal(m, m0)
+  expect_true(all(m$nsnp == 77))
+
+  expect_warning(
+    w <- mr_weighted_mode(d$beta.exposure, d$beta.outcome, d$se.exposure, d$se.outcome),
+    "Removing 2 SNP\\(s\\)"
+  )
+  expect_true(is.finite(w$b))
+  expect_equal(w$nsnp, 77)
+})
+
+test_that("mr_mode() returns NA when fewer than 3 SNPs have non-zero beta.exposure (issue #570)", {
+  d <- dat[1:4, ]
+  d$beta.exposure[1:2] <- 0
+  warnings <- character()
+  w <- withCallingHandlers(
+    mr_simple_mode(d$beta.exposure, d$beta.outcome, d$se.exposure, d$se.outcome),
+    warning = function(cnd) {
+      warnings <<- c(warnings, conditionMessage(cnd))
+      invokeRestart("muffleWarning")
+    }
+  )
+  expect_equal(w, list(b = NA, se = NA, pval = NA, nsnp = NA))
+  expect_match(warnings[1], "Removing 2 SNP\\(s\\)")
+  expect_match(warnings[2], "Need at least 3 SNPs")
+})

tests/testthat/test_mvmr_nsnp.R

@@ -0,0 +1,53 @@
+make_mvdat <- function() {
+  set.seed(1234)
+  nsnp <- 40
+  exposure_beta <- matrix(stats::rnorm(nsnp * 2, sd = 0.1), ncol = 2)
+  colnames(exposure_beta) <- c("exp1", "exp2")
+  # 30 SNPs pass 5e-8 for exp1 and 10 for exp2
+  exposure_pval <- cbind(
+    exp1 = c(rep(1e-10, 30), rep(0.01, 10)),
+    exp2 = c(rep(0.01, 30), rep(1e-10, 10))
+  )
+  outcome_beta <- as.vector(exposure_beta %*% c(0.5, -0.3)) + stats::rnorm(nsnp, sd = 0.01)
+  list(
+    exposure_beta = exposure_beta,
+    exposure_pval = exposure_pval,
+    exposure_se = matrix(0.01, nsnp, 2, dimnames = list(NULL, c("exp1", "exp2"))),
+    outcome_beta = outcome_beta,
+    outcome_se = rep(0.01, nsnp),
+    outcome_pval = rep(0.5, nsnp),
+    expname = data.frame(id.exposure = c("exp1", "exp2"), exposure = c("Exposure 1", "Exposure 2")),
+    outname = data.frame(id.outcome = "out", outcome = "Outcome")
+  )
+}
+
+test_that("mv_multiple() nsnp is all SNPs when instrument_specific = FALSE (issue #78)", {
+  mvdat <- make_mvdat()
+  res <- mv_multiple(mvdat)$result
+  expect_equal(res$nsnp, c(40, 40))
+  # pval_threshold does not affect the estimates when instrument_specific = FALSE
+  res1 <- mv_multiple(mvdat, pval_threshold = 1)$result
+  expect_equal(res, res1)
+  res_is <- mv_multiple(mvdat, instrument_specific = TRUE)$result
+  expect_equal(res_is$nsnp, c(30, 10))
+})
+
+test_that("mv_residual() nsnp is all SNPs when instrument_specific = FALSE (issue #78)", {
+  mvdat <- make_mvdat()
+  res <- mv_residual(mvdat)$result
+  expect_equal(res$nsnp, c(40, 40))
+  res1 <- mv_residual(mvdat, pval_threshold = 1)$result
+  expect_equal(res, res1)
+  res_is <- mv_residual(mvdat, instrument_specific = TRUE)$result
+  expect_equal(res_is$nsnp, c(30, 10))
+})
+
+test_that("mv_residual() estimates are not set to NA when few SNPs pass the threshold and instrument_specific = FALSE (issue #78)", {
+  mvdat <- make_mvdat()
+  # Only 1 SNP passes the threshold for exp2
+  mvdat$exposure_pval[, "exp2"] <- c(1e-10, rep(0.01, 39))
+  res <- mv_residual(mvdat)$result
+  expect_false(anyNA(res$b))
+  res_is <- mv_residual(mvdat, instrument_specific = TRUE)$result
+  expect_true(is.na(res_is$b[res_is$id.exposure == "exp2"]))
+})