From 55d176fca01072eb6890ced32527bc990137475e Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Mon, 19 Oct 2026 15:33:14 +0100 Subject: [PATCH 01/25] Fix URL Closes #205 --- vignettes/gwas2020.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/gwas2020.Rmd b/vignettes/gwas2020.Rmd index 8fcd7fbbf..b6c26b3a4 100644 --- a/vignettes/gwas2020.Rmd +++ b/vignettes/gwas2020.Rmd @@ -27,7 +27,7 @@ Another change is that the R package that managed the authentication has updated We conducted a large GWAS analysis using a pipeline that systematically analysed every [PHESANT](https://pubmed.ncbi.nlm.nih.gov/29040602/) phenotype in UK Biobank. There were previously ~20k traits with complete GWAS data, but a majority of these were binary traits based on very few numbers of cases. We have now filtered out unreliable datasets, there are 2514 traits remaining, with any binary traits removed that had fewer than 1000 cases. Another issue is the combination of small numbers of cases and allele frequency - here minor allele count (MAC) for a particular association could be very small which would lead to high false positives when using Bolt-LMM. The remaining traits have been filtered to only retain associations where the MAC > 90. -Document detailing this investigation here: https://htmlpreview.github.io/?https://raw.githubusercontent.com/MRCIEU/ukbb-gwas-analysis/master/docs/ldsc_clumped_analysis.html?token=AAOV6TBQXEXEPT7SUXXLWMC6DWP3O +Document detailing this investigation here: ### All data is now harmonised From a24632ac9054e026806f8d4165351ec886cad4fa Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:02:29 +0100 Subject: [PATCH 02/25] Bump version --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 9a1dfe23e..5777234f5 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.7.12 +Version: 0.7.13 Authors@R: c( person("Gibran", "Hemani", , "g.hemani@bristol.ac.uk", role = c("aut", "cre"), comment = c(ORCID = "0000-0003-0920-1055")), From fa0e57194c098c307470980cfc2f66926c3d9b39 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:02:38 +0100 Subject: [PATCH 03/25] Update NEWS.md --- NEWS.md | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/NEWS.md b/NEWS.md index 4e2e5e40c..56ffc7302 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,10 @@ +# TwoSampleMR v0.7.13 + +(Release date 2026-10-##) + +* 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) + # TwoSampleMR v0.7.12 (Release date 2026-10-08) From 0ea3984c83454653d399a3cf5700ac6179c1f32c Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:03:23 +0100 Subject: [PATCH 04/25] Fix outcome extraction object and clarify steps in MVMR section of Perform MR vignette Closes #340 --- vignettes/perform_mr.Rmd | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/vignettes/perform_mr.Rmd b/vignettes/perform_mr.Rmd index c4506e00d..44a4f1454 100644 --- a/vignettes/perform_mr.Rmd +++ b/vignettes/perform_mr.Rmd @@ -648,8 +648,11 @@ When SNPs instrument multiple potential exposures, for example in the case of di 2. Combine these into a set of all instruments 3. Clump these to avoid the possibility that e.g. a variant for exposure 1 is in LD with a variant for exposure 2 4. Re-extract all the final clumped SNPs from (3) from all of the exposures -5. Harmonise them all to be on the same effect allele -6. Use the multivariable MR method against these harmonised data +5. Extract the final clumped SNPs from (3) from the outcome +6. Harmonise them all to be on the same effect allele +7. Use the multivariable MR method against these harmonised data + +Steps 1 to 4 are performed by `mv_extract_exposures()` (or `mv_extract_exposures_local()` for your own summary data), step 5 by `extract_outcome_data()`, step 6 by `mv_harmonise_data()`, and step 7 by `mv_multiple()` or `mv_residual()`. Example - The GWAS IDs for HDL, LDL and total cholesterol are `ieu-a-299`, `ieu-a-300`, and `ieu-a-302`. The GWAS ID for coronary heart disease (CHD) is `ieu-a-7`. In this example we will estimate the multivariable effects of HDL, LDL and total cholesterol on CHD. @@ -667,7 +670,7 @@ mv_exposure_dat <- mv_extract_exposures(id_exposure) Next, also extract those SNPs from the outcome. ```{r eval=FALSE} -mv_outcome_dat <- extract_outcome_data(exposure_dat$SNP, id_outcome) +mv_outcome_dat <- extract_outcome_data(mv_exposure_dat$SNP, id_outcome) ``` Once the data has been obtained, harmonise so that all are on the same reference allele. From 4d47769cd25058aa211b16473e7f4e0d16d03984 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:07:36 +0100 Subject: [PATCH 05/25] Fix harmonise_ld_dat() when LD matrix SNPs are in a different order to the summary data Closes #348 --- R/other_formats.R | 10 +++++++--- tests/testthat/test_otherformats.R | 31 ++++++++++++++++++++++++++++++ 2 files changed, 38 insertions(+), 3 deletions(-) diff --git a/R/other_formats.R b/R/other_formats.R index 1c3a21222..9c60f23df 100644 --- a/R/other_formats.R +++ b/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)) diff --git a/tests/testthat/test_otherformats.R b/tests/testthat/test_otherformats.R index 63d6c2db6..dff9aa799 100644 --- a/tests/testthat/test_otherformats.R +++ b/tests/testthat/test_otherformats.R @@ -29,6 +29,37 @@ test_that("MRInput with cor", { }) +test_that("harmonise_ld_dat handles LD matrix in a different SNP order", { + x <- data.frame( + SNP = c("rs1", "rs2", "rs3"), + effect_allele.exposure = c("A", "C", "G"), + other_allele.exposure = c("G", "T", "A"), + beta.exposure = c(0.1, 0.2, 0.3), + beta.outcome = c(0.01, 0.02, 0.03) + ) + # LD matrix in a different order to x, with rs1 on the opposite effect allele + ldnames <- c("rs3_G_A", "rs1_G_A", "rs2_C_T") + ld <- matrix( + c(1, 0.13, 0.23, + 0.13, 1, 0.12, + 0.23, 0.12, 1), + nrow = 3, + dimnames = list(ldnames, ldnames) + ) + out <- harmonise_ld_dat(x, ld) + expect_equal(out$x$SNP, c("rs1", "rs2", "rs3")) + expect_equal(rownames(out$ld), c("rs1", "rs2", "rs3")) + expect_equal(colnames(out$ld), c("rs1", "rs2", "rs3")) + expect_equal(out$ld["rs1", "rs2"], 0.12) + expect_equal(out$ld["rs1", "rs3"], 0.13) + expect_equal(out$ld["rs2", "rs3"], 0.23) + expect_equal(out$x$beta.exposure, c(-0.1, 0.2, 0.3)) + expect_equal(out$x$beta.outcome, c(-0.01, 0.02, 0.03)) + expect_equal(out$x$effect_allele.exposure, c("G", "C", "G")) + expect_equal(out$x$other_allele.exposure, c("A", "T", "A")) +}) + + test_that("mrpresso", { w <- expect_warning(run_mr_presso(dat)) expect_true(length(w) == 1) From 9afb53e11d57cc1859a83551eb28e457f985ea5f Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:07:38 +0100 Subject: [PATCH 06/25] Update NEWS.md --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index 56ffc7302..95afc7477 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,7 @@ * 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) # TwoSampleMR v0.7.12 From b49b13336cb234b5d32d53b72862b8c8e45ab026 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:11:10 +0100 Subject: [PATCH 07/25] Stop mr_report() rounding small p-values to 0 Closes #42 --- R/knit.R | 14 ++++++++++++++ tests/testthat/test_knit.R | 15 +++++++++++++++ 2 files changed, 29 insertions(+) create mode 100644 tests/testthat/test_knit.R diff --git a/R/knit.R b/R/knit.R index 1c5b724ec..38ecdec96 100644 --- a/R/knit.R +++ b/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)) } diff --git a/tests/testthat/test_knit.R b/tests/testthat/test_knit.R new file mode 100644 index 000000000..0b619f7ba --- /dev/null +++ b/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"))) +}) From f5f40a060cde29057384e3d116656997f516944e Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:11:15 +0100 Subject: [PATCH 08/25] Update NEWS.md --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index 95afc7477..c5a5c9861 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,7 @@ * 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) # TwoSampleMR v0.7.12 From b916f672d2741042c0753332428f8c540c8a5651 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:18:03 +0100 Subject: [PATCH 09/25] Bump roxygen2 --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 5777234f5..62a173174 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -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 From a3d7f1b182a4296147a51fa38c9de2eede7ffe25 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:18:40 +0100 Subject: [PATCH 10/25] Remove SNPs with eaf.exposure of 0 or 1 in run_mrmix() Closes #169 --- R/other_formats.R | 19 +++++++++++++++++ tests/testthat/test_otherformats.R | 33 ++++++++++++++++++++++++++++++ 2 files changed, 52 insertions(+) diff --git a/R/other_formats.R b/R/other_formats.R index 9c60f23df..967ff0efa 100644 --- a/R/other_formats.R +++ b/R/other_formats.R @@ -255,6 +255,7 @@ mr_ivw_radial <- function(b_exp, b_out, se_exp, se_out, parameters = default_par #' See 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 @@ -287,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, diff --git a/tests/testthat/test_otherformats.R b/tests/testthat/test_otherformats.R index dff9aa799..abe666217 100644 --- a/tests/testthat/test_otherformats.R +++ b/tests/testthat/test_otherformats.R @@ -60,6 +60,39 @@ test_that("harmonise_ld_dat handles LD matrix in a different SNP order", { }) +test_that("run_mrmix removes SNPs with eaf.exposure of 0 or 1", { + skip_if_not_installed("MRMix") + d <- dat + d$eaf.exposure[1] <- 1 + d$eaf.exposure[2] <- 0 + expect_warning( + expect_message(w <- run_mrmix(d)), + "Removing 2 SNP\\(s\\)" + ) + expect_message(w0 <- run_mrmix(dat[-(1:2), ])) + expect_equal(w, w0) +}) + +test_that("run_mrmix returns NULL when no SNPs have a valid eaf.exposure", { + skip_if_not_installed("MRMix") + d <- dat + d$eaf.exposure <- 1 + warnings <- character() + expect_message( + w <- withCallingHandlers( + run_mrmix(d), + warning = function(cnd) { + warnings <<- c(warnings, conditionMessage(cnd)) + invokeRestart("muffleWarning") + } + ) + ) + expect_match(warnings[1], "Removing 79 SNP\\(s\\)") + expect_match(warnings[2], "No SNPs remaining") + expect_null(w[[1]]) +}) + + test_that("mrpresso", { w <- expect_warning(run_mr_presso(dat)) expect_true(length(w) == 1) From 981c67e11b66e5e9ff00592abf807df9e8bf3a07 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:18:47 +0100 Subject: [PATCH 11/25] Document --- man/run_mrmix.Rd | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/man/run_mrmix.Rd b/man/run_mrmix.Rd index 800e9814a..f0ed463e2 100644 --- a/man/run_mrmix.Rd +++ b/man/run_mrmix.Rd @@ -7,7 +7,8 @@ run_mrmix(dat) } \arguments{ -\item{dat}{Output from \code{\link[=harmonise_data]{harmonise_data()}}. Ensures that no eaf.exposure values are missing.} +\item{dat}{Output from \code{\link[=harmonise_data]{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.} } \value{ List of results, with one list item for every exposure/outcome pair in dat object From 0fa07201d41eb29b86c50669f17cea9ebdcaa909 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:18:50 +0100 Subject: [PATCH 12/25] Update NEWS.md --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index c5a5c9861..2ebc3f33e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -6,6 +6,7 @@ * 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) # TwoSampleMR v0.7.12 From 3af456696eaad780afcc7b145163d772429ad2cc Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:23:52 +0100 Subject: [PATCH 13/25] Remove SNPs with beta.exposure of 0 in mode based estimators Closes #570 --- R/mr_mode.R | 17 ++++++++++++++++- tests/testthat/test_mr_mode.R | 34 ++++++++++++++++++++++++++++++++++ 2 files changed, 50 insertions(+), 1 deletion(-) diff --git a/R/mr_mode.R b/R/mr_mode.R index 4085d0d25..06461c805 100644 --- a/R/mr_mode.R +++ b/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 diff --git a/tests/testthat/test_mr_mode.R b/tests/testthat/test_mr_mode.R index 2483b7ac7..1e0d4ae54 100644 --- a/tests/testthat/test_mr_mode.R +++ b/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") +}) From 1da60513ea2958cf97a09dde9a7f16ee95f12eb7 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:23:56 +0100 Subject: [PATCH 14/25] Update NEWS.md --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index 2ebc3f33e..1c5ea86ee 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,6 +7,7 @@ * 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) # TwoSampleMR v0.7.12 From 72312e88839d212d92c856cb3db4e56f68dcf68e Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:26:24 +0100 Subject: [PATCH 15/25] Fix nsnp in mv_multiple() and mv_residual() when instrument_specific = FALSE Closes #78 --- R/multivariable_mr.R | 22 ++++++++++---- tests/testthat/test_mvmr_nsnp.R | 53 +++++++++++++++++++++++++++++++++ 2 files changed, 69 insertions(+), 6 deletions(-) create mode 100644 tests/testthat/test_mvmr_nsnp.R diff --git a/R/multivariable_mr.R b/R/multivariable_mr.R index 25f13387a..32298b77e 100644 --- a/R/multivariable_mr.R +++ b/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 diff --git a/tests/testthat/test_mvmr_nsnp.R b/tests/testthat/test_mvmr_nsnp.R new file mode 100644 index 000000000..976f89955 --- /dev/null +++ b/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"])) +}) From eb5ebc1d4f06804190d6d97ff630cdb35dba556c Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:26:30 +0100 Subject: [PATCH 16/25] Document --- man/mv_multiple.Rd | 2 +- man/mv_residual.Rd | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/man/mv_multiple.Rd b/man/mv_multiple.Rd index e705edbb5..c73984254 100644 --- a/man/mv_multiple.Rd +++ b/man/mv_multiple.Rd @@ -19,7 +19,7 @@ mv_multiple( \item{instrument_specific}{Should the estimate for each exposure be obtained by using all instruments from all exposures (\code{FALSE}, default) or by using only the instruments specific to each exposure (\code{TRUE}).} -\item{pval_threshold}{P-value threshold to include instruments. The default is \code{5e-8}.} +\item{pval_threshold}{P-value threshold to include instruments. Only used when \code{instrument_specific = TRUE}. The default is \code{5e-8}.} \item{plots}{Create plots? The default is \code{FALSE}.} } diff --git a/man/mv_residual.Rd b/man/mv_residual.Rd index ff4230c89..c4a913c47 100644 --- a/man/mv_residual.Rd +++ b/man/mv_residual.Rd @@ -19,7 +19,7 @@ mv_residual( \item{instrument_specific}{Should the estimate for each exposure be obtained by using all instruments from all exposures (\code{FALSE}, default) or by using only the instruments specific to each exposure (\code{TRUE}).} -\item{pval_threshold}{P-value threshold to include instruments. The default is \code{5e-8}.} +\item{pval_threshold}{P-value threshold to include instruments. Only used when \code{instrument_specific = TRUE}. The default is \code{5e-8}.} \item{plots}{Create plots? The default is \code{FALSE}.} } From d090851d6124eed9d11597d8e77b0c845de73fa4 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:26:34 +0100 Subject: [PATCH 17/25] Update NEWS.md --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index 1c5ea86ee..fbdd80de0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,6 +8,7 @@ * 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) # TwoSampleMR v0.7.12 From 09e160ae72b0e58024be6bfff60b66df23419198 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:37:00 +0100 Subject: [PATCH 18/25] Document the effective_n() formula and how it differs from METAL Closes #433 --- R/add_rsq.r | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/R/add_rsq.r b/R/add_rsq.r index 818b5eaef..01e035226 100644 --- a/R/add_rsq.r +++ b/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 +#' 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) . +#' +#' 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 From de660ee406b65735a2d25d4fb966b3b3c56269c1 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:37:07 +0100 Subject: [PATCH 19/25] Document --- man/effective_n.Rd | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/man/effective_n.Rd b/man/effective_n.Rd index 598cdaafd..e1b589962 100644 --- a/man/effective_n.Rd +++ b/man/effective_n.Rd @@ -15,5 +15,16 @@ effective_n(ncase, ncontrol) Vector of effective sample size } \description{ -Taken from \url{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) \url{https://www.nature.com/articles/nprot.2014.071}. +} +\details{ +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 \code{\link[=add_rsq]{add_rsq()}} to obtain the sample sizes for binary traits (those with units of \code{"log odds"}) which are then used by \code{\link[=directionality_test]{directionality_test()}} and \code{\link[=steiger_filtering]{steiger_filtering()}}. } From 86770884ac9131ed9c549c14fb9545b51cf9c289 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:37:11 +0100 Subject: [PATCH 20/25] Update NEWS.md --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index fbdd80de0..aeef0d400 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,6 +9,7 @@ * `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 From f2a2b651abaf87aa377752a4237fa94793aaae7a Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:49:18 +0100 Subject: [PATCH 21/25] Update NEWS.md --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index aeef0d400..25d7a9b33 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,6 @@ # TwoSampleMR v0.7.13 -(Release date 2026-10-##) +(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) From fbf4824296631fad76e91691fe43685483680a5f Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:59:42 +0100 Subject: [PATCH 22/25] Add urlcheck recipe --- justfile | 2 ++ 1 file changed, 2 insertions(+) diff --git a/justfile b/justfile index 9ebf294e6..737733832 100644 --- a/justfile +++ b/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()" From 18dc492ee7a4b66fbb895724de02bb6ddb045d03 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 11:59:55 +0100 Subject: [PATCH 23/25] Use doi --- inst/CITATION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inst/CITATION b/inst/CITATION index 30a8f75fb..2024806e6 100644 --- a/inst/CITATION +++ b/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( From ec306a9153a41fabb4034c8a5223d4adc019f458 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 12:00:06 +0100 Subject: [PATCH 24/25] Quote to make old URL not a URL --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 25d7a9b33..50317ed11 100644 --- a/NEWS.md +++ b/NEWS.md @@ -82,7 +82,7 @@ (Release date 2026-05-04) -* Amend URLs to +* Amend `https://gwas.mrcieu.ac.uk` URLs to * 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 From ae6f7d4b90973da7104d8f25636bb22a7dceb6d7 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Tue, 20 Oct 2026 12:00:13 +0100 Subject: [PATCH 25/25] Use DOI --- vignettes/gwas2020.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/gwas2020.Rmd b/vignettes/gwas2020.Rmd index b6c26b3a4..acbf7f28d 100644 --- a/vignettes/gwas2020.Rmd +++ b/vignettes/gwas2020.Rmd @@ -25,7 +25,7 @@ Another change is that the R package that managed the authentication has updated ### UKBiobank data has been curated -We conducted a large GWAS analysis using a pipeline that systematically analysed every [PHESANT](https://pubmed.ncbi.nlm.nih.gov/29040602/) phenotype in UK Biobank. There were previously ~20k traits with complete GWAS data, but a majority of these were binary traits based on very few numbers of cases. We have now filtered out unreliable datasets, there are 2514 traits remaining, with any binary traits removed that had fewer than 1000 cases. Another issue is the combination of small numbers of cases and allele frequency - here minor allele count (MAC) for a particular association could be very small which would lead to high false positives when using Bolt-LMM. The remaining traits have been filtered to only retain associations where the MAC > 90. +We conducted a large GWAS analysis using a pipeline that systematically analysed every [PHESANT](https://doi.org/10.1093/ije/dyx204) phenotype in UK Biobank. There were previously ~20k traits with complete GWAS data, but a majority of these were binary traits based on very few numbers of cases. We have now filtered out unreliable datasets, there are 2514 traits remaining, with any binary traits removed that had fewer than 1000 cases. Another issue is the combination of small numbers of cases and allele frequency - here minor allele count (MAC) for a particular association could be very small which would lead to high false positives when using Bolt-LMM. The remaining traits have been filtered to only retain associations where the MAC > 90. Document detailing this investigation here: