diff --git a/DESCRIPTION b/DESCRIPTION
index 9a1dfe23e..62a173174 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")),
@@ -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
diff --git a/NEWS.md b/NEWS.md
index 4e2e5e40c..50317ed11 100644
--- a/NEWS.md
+++ b/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
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 <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
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/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/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/R/other_formats.R b/R/other_formats.R
index 1c3a21222..967ff0efa 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))
@@ -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,
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(
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()"
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()}}.
 }
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}.}
 }
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
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")))
+})
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")
+})
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"]))
+})
diff --git a/tests/testthat/test_otherformats.R b/tests/testthat/test_otherformats.R
index 63d6c2db6..abe666217 100644
--- a/tests/testthat/test_otherformats.R
+++ b/tests/testthat/test_otherformats.R
@@ -29,6 +29,70 @@ 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("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)
diff --git a/vignettes/gwas2020.Rmd b/vignettes/gwas2020.Rmd
index 8fcd7fbbf..acbf7f28d 100644
--- a/vignettes/gwas2020.Rmd
+++ b/vignettes/gwas2020.Rmd
@@ -25,9 +25,9 @@ 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: 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: <https://htmlpreview.github.io/?https://raw.githubusercontent.com/MRCIEU/ukbb-gwas-analysis/master/docs/ldsc_clumped_analysis.html>
 
 ### All data is now harmonised
 
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.