Merge pull request #679 from remlapmot/devel-2026-04-30

remlapmot · web-flow · commit 48bed40d3c95 · 2026-05-04T11:07:54.000+01:00
TwoSampleMR 0.7.5
diff --git a/.Rbuildignore b/.Rbuildignore
@@ -21,3 +21,5 @@ cache$
 ^inst/sandpit$
 ^air\.toml$
 ^vignettes/rf\.rdata$
+^\.positai$
+^\.claude$
diff --git a/.github/workflows/check-full.yaml b/.github/workflows/check-full.yaml
@@ -38,6 +38,7 @@ jobs:
           - {os: ubuntu-latest,   r: '4.3.2'}
           - {os: ubuntu-latest,   r: 'oldrel-3'}
           - {os: ubuntu-latest,   r: 'oldrel-4'}
+          - {os: ubuntu-latest,   r: 'oldrel-5'}  
           - {os: macos-14, r: 'release'}
           # - {os: ubuntu-24.04-arm, r: 'release', rspm: 'no' }
 
diff --git a/.gitignore b/.gitignore
@@ -26,3 +26,4 @@ vignettes/mr_report.md
 vignettes/figure/*.png
 vignettes/*.html
 rf.rdata
+.positai
diff --git 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.4
+Version: 0.7.5
 Authors@R: c(
     person("Gibran", "Hemani", , "g.hemani@bristol.ac.uk", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0003-0920-1055")),
@@ -20,7 +20,7 @@ Authors@R: c(
   )
 Description: A package for performing Mendelian randomization using GWAS
     summary data. It uses the IEU OpenGWAS database
-    <https://gwas.mrcieu.ac.uk/> to automatically obtain data, and a wide
+    <https://opengwas.io> to automatically obtain data, and a wide
     range of methods to run the analysis.
 License: MIT + file LICENSE
 URL: https://github.com/MRCIEU/TwoSampleMR,
@@ -68,6 +68,6 @@ Remotes:
     MRPRESSO=rondolab/MR-PRESSO,
     qingyuanzhao/mr.raps,
     WSpiller/RadialMR
+Config/roxygen2/markdown: TRUE
+Config/roxygen2/version: 8.0.0
 Encoding: UTF-8
-Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.3
diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,15 @@
+# TwoSampleMR v0.7.5
+
+(Release date 2026-05-04)
+
+* 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
+* Pre-generate jackknife resample indices in `mr_rucker_jackknife_internal` and switch to `lapply`, matching the bootstrap function pattern
+* Vectorise `get_population_allele_frequency` by inlining the quadratic-root logic from `contingency()` and picking the valid root with vectorised constraint checks
+* Skip `mr_wald_ratio` inside `mr()` when more than one SNP is provided and other methods are also requested, suppressing the multi-SNP warning for default method runs while preserving the warning when `mr_wald_ratio` is requested on its own
+
 # TwoSampleMR v0.7.4
 
 (Release date 2026-04-02)
diff --git a/R/add_rsq.r b/R/add_rsq.r
@@ -347,13 +347,35 @@ allele_frequency <- function(g) {
 get_population_allele_frequency <- function(af, prop, odds_ratio, prevalence) {
   stopifnot(length(af) == length(odds_ratio))
   stopifnot(length(prop) == length(odds_ratio))
-  for (i in seq_along(odds_ratio)) {
-    co <- contingency(af[i], prop[i], odds_ratio[i])
-    af_controls <- co[1, 2] / (co[1, 2] + co[2, 2])
-    af_cases <- co[1, 1] / (co[1, 1] + co[2, 1])
-    af[i] <- af_controls * (1 - prevalence[i]) + af_cases * prevalence[i]
+  eps <- 1e-15
+  a <- odds_ratio - 1
+  b <- (af + prop) * (1 - odds_ratio) - 1
+  c_ <- odds_ratio * af * prop
+
+  z <- numeric(length(odds_ratio))
+  linear <- abs(a) < eps
+  z[linear] <- -c_[linear] / b[linear]
+
+  quad <- !linear
+  if (any(quad)) {
+    d <- pmax(0, b[quad]^2 - 4 * a[quad] * c_[quad])
+    sqrt_d <- sqrt(d)
+    two_a <- 2 * a[quad]
+    z_pos <- (-b[quad] + sqrt_d) / two_a
+    z_neg <- (-b[quad] - sqrt_d) / two_a
+    af_q <- af[quad]
+    prop_q <- prop[quad]
+    tol <- -1e-7
+    valid_pos <- z_pos >= tol &
+      (prop_q - z_pos) >= tol &
+      (af_q - z_pos) >= tol &
+      (1 + z_pos - af_q - prop_q) >= tol
+    z[quad] <- ifelse(valid_pos, z_pos, z_neg)
   }
-  return(af)
+
+  af_controls <- (af - z) / (1 - prop)
+  af_cases <- z / prop
+  return(af_controls * (1 - prevalence) + af_cases * prevalence)
 }
 
 
diff --git a/R/leaveoneout.R b/R/leaveoneout.R
@@ -44,7 +44,7 @@ mr_leaveoneout <- function(dat, parameters = default_parameters(), method = mr_i
       return(d)
     }
     if (nsnp > 2) {
-      l <- lapply(1:nsnp, function(i) {
+      l <- lapply(seq_len(nsnp), function(i) {
         with(
           x,
           method(beta.exposure[-i], beta.outcome[-i], se.exposure[-i], se.outcome[-i], parameters)
diff --git a/R/mr.R b/R/mr.R
@@ -72,15 +72,21 @@ mr <- function(
     } else {
       message("Analysing '", exp_id, "' on '", out_id, "'")
     }
-    res <- lapply(method_list, function(meth) {
+    keep <- if (nrow(x) > 1 && length(method_list) > 1) {
+      method_list != "mr_wald_ratio"
+    } else {
+      rep(TRUE, length(method_list))
+    }
+    methods_to_run <- method_list[keep]
+    res <- lapply(methods_to_run, function(meth) {
       get(meth)(x$beta.exposure, x$beta.outcome, x$se.exposure, x$se.outcome, parameters)
     })
     mr_tab <- data.frame(
       id.exposure = exp_id,
       id.outcome = out_id,
       outcome = x$outcome[1],
       exposure = x$exposure[1],
-      method = method_names,
+      method = method_names[keep],
       nsnp = vapply(res, function(x) x$nsnp, numeric(1)),
       b = vapply(res, function(x) x$b, numeric(1)),
       se = vapply(res, function(x) x$se, numeric(1)),
diff --git a/R/mr_mode.R b/R/mr_mode.R
@@ -78,7 +78,7 @@ mr_mode <- function(dat, parameters = default_parameters(), mode_method = "all")
     #Set up a matrix to store the results from each bootstrap iteration
     beta.boot <- matrix(nrow = nboot, ncol = length(beta_Mode.in))
 
-    for (i in 1:nboot) {
+    for (i in seq_len(nboot)) {
       BetaIV.boot <- BetaIV.boot_mat[i, ]
       BetaIV.boot_NOME <- BetaIV.boot_NOME_mat[i, ]
 
diff --git a/R/multivariable_mr.R b/R/multivariable_mr.R
@@ -498,15 +498,15 @@ mv_residual <- function(
   pval.exposure <- mvdat$exposure_pval
 
   nexp <- ncol(beta.exposure)
-  effs <- array(1:nexp)
-  se <- array(1:nexp)
-  pval <- array(1:nexp)
-  nsnp <- array(1:nexp)
+  effs <- numeric(nexp)
+  se <- numeric(nexp)
+  pval <- numeric(nexp)
+  nsnp <- integer(nexp)
   marginal_outcome <- matrix(0, nrow(beta.exposure), ncol(beta.exposure))
   p <- list()
   nom <- colnames(beta.exposure)
   nom2 <- mvdat$expname$exposure[match(nom, mvdat$expname$id.exposure)]
-  for (i in 1:nexp) {
+  for (i in seq_len(nexp)) {
     # For this exposure, only keep SNPs that meet some p-value threshold
     index <- pval.exposure[, i] < pval_threshold
 
@@ -607,15 +607,15 @@ mv_multiple <- function(
   w <- 1 / mvdat$outcome_se^2
 
   nexp <- ncol(beta.exposure)
-  effs <- array(1:nexp)
-  se <- array(1:nexp)
-  pval <- array(1:nexp)
-  nsnp <- array(1:nexp)
+  effs <- numeric(nexp)
+  se <- numeric(nexp)
+  pval <- numeric(nexp)
+  nsnp <- integer(nexp)
   # marginal_outcome <- matrix(0, nrow(beta.exposure), ncol(beta.exposure))
   p <- list()
   nom <- colnames(beta.exposure)
   nom2 <- mvdat$expname$exposure[match(nom, mvdat$expname$id.exposure)]
-  for (i in 1:nexp) {
+  for (i in seq_len(nexp)) {
     # For this exposure, only keep SNPs that meet some p-value threshold
     index <- pval.exposure[, i] < pval_threshold
 
@@ -704,15 +704,15 @@ mv_basic <- function(mvdat, pval_threshold = 5e-8) {
   pval.exposure <- mvdat$exposure_pval
 
   nexp <- ncol(beta.exposure)
-  effs <- array(1:nexp)
-  se <- array(1:nexp)
-  pval <- array(1:nexp)
-  nsnp <- array(1:nexp)
+  effs <- numeric(nexp)
+  se <- numeric(nexp)
+  pval <- numeric(nexp)
+  nsnp <- integer(nexp)
   marginal_outcome <- matrix(0, nrow(beta.exposure), ncol(beta.exposure))
   p <- list()
   nom <- colnames(beta.exposure)
   nom2 <- mvdat$expname$exposure[match(nom, mvdat$expname$id.exposure)]
-  for (i in 1:nexp) {
+  for (i in seq_len(nexp)) {
     # For this exposure, only keep SNPs that meet some p-value threshold
     index <- pval.exposure[, i] < pval_threshold
 
@@ -772,15 +772,15 @@ mv_ivw <- function(mvdat, pval_threshold = 5e-8) {
   w <- 1 / mvdat$outcome_se^2
 
   nexp <- ncol(beta.exposure)
-  effs <- array(1:nexp)
-  se <- array(1:nexp)
-  pval <- array(1:nexp)
-  nsnp <- array(1:nexp)
+  effs <- numeric(nexp)
+  se <- numeric(nexp)
+  pval <- numeric(nexp)
+  nsnp <- integer(nexp)
   # marginal_outcome <- matrix(0, nrow(beta.exposure), ncol(beta.exposure))
   p <- list()
   nom <- colnames(beta.exposure)
   nom2 <- mvdat$expname$exposure[match(nom, mvdat$expname$id.exposure)]
-  for (i in 1:nexp) {
+  for (i in seq_len(nexp)) {
     # For this exposure, only keep SNPs that meet some p-value threshold
     index <- pval.exposure[, i] < pval_threshold
 
diff --git a/R/rucker.R b/R/rucker.R
@@ -37,7 +37,7 @@ PM <- function(y = y, s = s, Alpha = 0.1) {
   v <- 1 / s^2
   sum.v <- sum(v)
   typS <- sum(v * (k - 1)) / (sum.v^2 - sum(v^2))
-  for (j in 1:L) {
+  for (j in seq_len(L)) {
     tausq <- 0
     Fstat <- 1
     TAUsq <- NULL
@@ -302,7 +302,7 @@ mr_rucker_bootstrap <- function(dat, parameters = default_parameters()) {
 
   dat2 <- dat
   l <- list()
-  for (i in 1:nboot) {
+  for (i in seq_len(nboot)) {
     dat2$beta.exposure <- boot_exp[i, ]
     dat2$beta.outcome <- boot_out[i, ]
     l[[i]] <- mr_rucker(dat2, parameters)
@@ -451,13 +451,15 @@ mr_rucker_jackknife_internal <- function(dat, parameters = default_parameters())
       e_plot = NULL
     ))
   } else {
-    l <- list()
-    for (i in 1:nboot) {
-      # dat2$beta.exposure <- rnorm(nsnp, mean=dat$beta.exposure, sd=dat$se.exposure)
-      # dat2$beta.outcome <- rnorm(nsnp, mean=dat$beta.outcome, sd=dat$se.outcome)
-      dat2 <- dat[sample(seq_len(nrow(dat)), nrow(dat), replace = TRUE), ]
-      l[[i]] <- mr_rucker_internal(dat2, parameters)
-    }
+    n <- nrow(dat)
+    boot_idx <- matrix(
+      sample.int(n, n * nboot, replace = TRUE),
+      nrow = nboot,
+      ncol = n
+    )
+    l <- lapply(seq_len(nboot), function(i) {
+      mr_rucker_internal(dat[boot_idx[i, ], ], parameters)
+    })
 
     modsel <- data.table::rbindlist(
       lapply(l, function(x) x$selected),
@@ -505,7 +507,7 @@ mr_rucker_jackknife_internal <- function(dat, parameters = default_parameters())
 
     p1 <- ggplot2::ggplot(bootstrap, ggplot2::aes_string(x = "Q", y = "Qdash")) +
       ggplot2::geom_point(ggplot2::aes_string(colour = "model")) +
-      ggplot2::geom_point(data = subset(bootstrap, i == "Full")) +
+      ggplot2::geom_point(data = bootstrap[bootstrap$i == "Full", ]) +
       ggplot2::scale_colour_brewer(type = "qual") +
       ggplot2::xlim(0, max(bootstrap$Q, bootstrap$Qdash)) +
       ggplot2::ylim(0, max(bootstrap$Q, bootstrap$Qdash)) +
diff --git a/R/singlesnp.R b/R/singlesnp.R
@@ -59,7 +59,7 @@ mr_singlesnp <- function(
       )
       return(d)
     }
-    l <- lapply(1:nsnp, function(j) {
+    l <- lapply(seq_len(nsnp), function(j) {
       with(
         x,
         get(single_method)(
diff --git a/README.md b/README.md
@@ -16,7 +16,7 @@ badge](https://mrcieu.r-universe.dev/badges/TwoSampleMR)](https://mrcieu.r-unive
 <!-- badges: end -->
 
 A package for performing Mendelian randomization using GWAS summary
-data. It uses the [IEU OpenGWAS database](https://gwas.mrcieu.ac.uk/) to
+data. It uses the [IEU OpenGWAS database](https://opengwas.io) to
 obtain data automatically, and a wide range of methods to run the
 analysis.
 
diff --git a/index.md b/index.md
@@ -7,7 +7,7 @@
 [![TwoSampleMR status badge](https://mrcieu.r-universe.dev/badges/TwoSampleMR)](https://mrcieu.r-universe.dev/TwoSampleMR)
 <!-- badges: end -->
 
-A package for performing Mendelian randomization using GWAS summary data. It uses the [IEU OpenGWAS database](https://gwas.mrcieu.ac.uk/) to obtain data automatically, and a wide range of methods to run the analysis.
+A package for performing Mendelian randomization using GWAS summary data. It uses the [IEU OpenGWAS database](https://opengwas.io) to obtain data automatically, and a wide range of methods to run the analysis.
 
 ## January 2020 major update 
 
diff --git a/man/TwoSampleMR-package.Rd b/man/TwoSampleMR-package.Rd
diff --git a/man/mr_steiger.Rd b/man/mr_steiger.Rd
diff --git a/man/mr_steiger2.Rd b/man/mr_steiger2.Rd
diff --git a/man/steiger_sensitivity.Rd b/man/steiger_sensitivity.Rd
diff --git a/vignettes/exposure.Rmd b/vignettes/exposure.Rmd
@@ -281,7 +281,7 @@ cg25212131_exp_dat <-
 
 ### IEU OpenGWAS database
 
-The IEU OpenGWAS database contains the entire summary statistics for thousands of GWASs. You can browse them here: https://gwas.mrcieu.ac.uk/
+The IEU OpenGWAS database contains the entire summary statistics for thousands of GWASs. You can browse them here: https://opengwas.io
 
 You can use this database to define the instruments for a particular exposure. You can also use this database to obtain the effects for constructing polygenic risk scores using different p-value thresholds.
 
diff --git a/vignettes/gwas2020.Rmd b/vignettes/gwas2020.Rmd
@@ -53,7 +53,7 @@ We are using Elasticsearch and Neo4j on an Oracle Cloud Infrastructure to serve
 
 ### Browse available datasets online
 
-We have a new home for the GWAS summary data: <https://gwas.mrcieu.ac.uk/>.
+We have a new home for the GWAS summary data: <https://opengwas.io>.
 
 ### Chromosome and position
 
@@ -89,15 +89,15 @@ It is now possible to perform clumping, or create LD matrices, using your own lo
 
 ### Access the data directly
 
-Previously the data was only accessible through the database. Now the data can be downloaded in "GWAS VCF" format from here https://gwas.mrcieu.ac.uk/. (IEU members can access all the data on RDSF or bluecrystal4 directly). This means that if you want to perform very large or numerous operations, you can do it on HPC or locally in a more performant manner by using the data files directly. Please see the [gwasvcf R package](https://github.com/mrcieu/gwasvcf) on how to work with these data.
+Previously the data was only accessible through the database. Now the data can be downloaded in "GWAS VCF" format from here https://opengwas.io. (IEU members can access all the data on RDSF or bluecrystal4 directly). This means that if you want to perform very large or numerous operations, you can do it on HPC or locally in a more performant manner by using the data files directly. Please see the [gwasvcf R package](https://github.com/mrcieu/gwasvcf) on how to work with these data.
 
 ### Connect the data to different analytical tools
 
 Either the data in the database, or the GWAS VCF files, can be queried and the results translated into the formats for a bunch of different R packages for MR, colocalisation, fine mapping, etc. Have a look at the [gwasglue R package](https://github.com/mrcieu/gwasglue), to see what is available and how to do this. It's still under construction, but feel free to try it, make suggestions, and contribute code.
 
 ## Key links
 
-- The IEU OpenGWAS database: https://gwas.mrcieu.ac.uk
+- The IEU OpenGWAS database: https://opengwas.io
 - API to the IEU OpenGWAS database: https://api.opengwas.io/api/
 - ieugwasr package, for R access to the API: https://mrcieu.github.io/ieugwasr/
 - ieugwaspy package, for python access to the API: https://github.com/MRCIEU/ieugwaspy/
diff --git a/vignettes/introduction.Rmd b/vignettes/introduction.Rmd
@@ -82,7 +82,7 @@ install.packages('TwoSampleMR',
 The workflow for performing MR is as follows:
 
 1. Select instruments for the exposure (perform LD clumping if necessary)
-2. Extract the instruments from the [IEU GWAS database](https://gwas.mrcieu.ac.uk/) for the outcomes of interest
+2. Extract the instruments from the [IEU GWAS database](https://opengwas.io) for the outcomes of interest
 3. Harmonise the effect sizes for the instruments on the exposures and the outcomes to be each for the same reference allele
 4. Perform MR analysis, sensitivity analyses, create plots, compile reports
 
diff --git a/vignettes/outcome.Rmd b/vignettes/outcome.Rmd

Original file line number	Diff line number	Diff line change
`@@ -44,7 +44,7 @@ mr_leaveoneout <- function(dat, parameters = default_parameters(), method = mr_i`
`44`	`44`	`return(d)`
`45`	`45`	`}`
`46`	`46`	`if (nsnp > 2) {`
`47`		`- l <- lapply(1:nsnp, function(i) {`
	`47`	`+ l <- lapply(seq_len(nsnp), function(i) {`
`48`	`48`	`with(`
`49`	`49`	`x,`
`50`	`50`	`method(beta.exposure[-i], beta.outcome[-i], se.exposure[-i], se.outcome[-i], parameters)`