diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index 77ee09f2..07c7735e 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.github/workflows/test-coverage.yaml @@ -49,7 +49,7 @@ jobs: covr::to_cobertura(cov) shell: Rscript {0} - - uses: codecov/codecov-action@v6 + - uses: codecov/codecov-action@v7 with: # Fail if error if not on PR, or if on PR and token is given fail_ci_if_error: ${{ github.event_name != 'pull_request' || secrets.CODECOV_TOKEN }} diff --git a/DESCRIPTION b/DESCRIPTION index 000d8749..03200344 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.7 +Version: 0.7.8 Authors@R: c( person("Gibran", "Hemani", , "g.hemani@bristol.ac.uk", role = c("aut", "cre"), comment = c(ORCID = "0000-0003-0920-1055")), diff --git a/NEWS.md b/NEWS.md index 100dbb4c..ad040606 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,20 @@ +# TwoSampleMR v0.7.8 + +(Release date 2026-06-10) + +* Fixed a copy-paste bug in `ldsc_rg()` which passed `d1$l2` instead of `d2$l2` as the weight vector to the trait-2 `ldsc_h2_internal()` call +* Fixed `extract_split()` chunking which gave `nsplit = 0` for SNP lists smaller than half the `splitsize`; now uses `ceiling()` with contiguous chunk ids +* Fixed the penalised mode in `mr_mode()` for vector-valued `phi`: each SNP is now penalised against the weighted-mode estimate for the same `phi` rather than recycling the mode vector across SNPs (results unchanged for the default scalar `phi = 1`) +* Fixed row order restoration in `add_metadata()` to use `dat[order(dat[[order_col]]), ]` +* Qualified a bare `aes()` call as `ggplot2::aes()` in `test_r_from_pn()` +* Used `isTRUE(all.equal())` instead of `all.equal() == TRUE` inside `ifelse()` for null-line detection in the forest plot functions, which could otherwise return a multi-element vector +* Simplified the group `Index` assignment in `forest_plot_1_to_many()` using `Letters[match()]`, which is also robust to non-contiguous group rows +* Collapsed a duplicated `||` condition in `extract_outcome_data_internal()` +* Removed dead `TAUsq` accumulation in `PM()` in `rucker.R` +* Fixed an error message typo in `combine_data()` and switched its `sapply()` to `vapply()` +* Replaced superseded `tidyr::gather()` with `tidyr::pivot_longer()` in `test_r_from_pn()` and dropped a redundant `requireNamespace("tidyr")` guard +* Lint package with Jarl + # TwoSampleMR v0.7.7 (Release date 2026-06-07) diff --git a/R/add_metadata.r b/R/add_metadata.r index c50acf5b..c9b04923 100644 --- a/R/add_metadata.r +++ b/R/add_metadata.r @@ -75,7 +75,7 @@ add_metadata <- function(dat, cols = c("sample_size", "ncase", "ncontrol", "unit names(dat)[names(dat) == "unit.exposure"] <- "units.exposure" names(dat)[names(dat) == "unit.outcome"] <- "units.outcome" - dat <- dat[dat[[order_col]], ] + dat <- dat[order(dat[[order_col]]), ] dat <- dat[, !names(dat) %in% order_col] # dat <- fix_ukb_d(dat) return(dat) diff --git a/R/add_rsq.r b/R/add_rsq.r index ec333821..818b5eae 100644 --- a/R/add_rsq.r +++ b/R/add_rsq.r @@ -73,7 +73,7 @@ add_rsq_one <- function(dat, what = "exposure") { message("Try adding metadata with add_metadata()") } } else if ( - all(grepl("SD", dat[[paste0("units.", what)]])) && all(!is.na(dat[[paste0("eaf.", what)]])) + all(grepl("SD", dat[[paste0("units.", what)]])) && !anyNA(dat[[paste0("eaf.", what)]]) ) { dat[[paste0("rsq.", what)]] <- NA dat[[paste0("rsq.", what)]] <- 2 * @@ -115,13 +115,6 @@ get_r_from_pn_less_accurate <- function(p, n) { } test_r_from_pn <- function() { - if (!requireNamespace("tidyr", quietly = TRUE)) { - stop( - "Package \"tidyr\" must be installed to use this function.", - call. = FALSE - ) - } - param <- expand.grid( n = c(10, 100, 1000, 10000, 100000), rsq = 10^seq(-4, -0.5, length.out = 30) @@ -137,11 +130,16 @@ test_r_from_pn <- function() { param$rsq2[i] <- (get_r_from_pn(param$pval[i], param$n[i]))^2 } - param <- tidyr::gather(param, key = out, value = value, rsq1, rsq2) + param <- tidyr::pivot_longer( + param, + cols = c("rsq1", "rsq2"), + names_to = "out", + values_to = "value" + ) p <- ggplot2::ggplot(param, ggplot2::aes(x = rsq_emp, value)) + ggplot2::geom_abline(slope = 1, linetype = "dotted") + - ggplot2::geom_line(aes(colour = out)) + + ggplot2::geom_line(ggplot2::aes(colour = out)) + ggplot2::facet_grid(. ~ n) + ggplot2::scale_x_log10() + ggplot2::scale_y_log10() diff --git a/R/forest_plot.R b/R/forest_plot.R index 2523d76d..676da589 100644 --- a/R/forest_plot.R +++ b/R/forest_plot.R @@ -323,7 +323,7 @@ mr_forest_plot_grouped <- } ## act on debug flag - if (debug == FALSE) { + if (!debug) { options(warn = -1) } @@ -372,7 +372,7 @@ mr_forest_plot_grouped <- ## draw and export the annotated forest plot grid.newpage() grid.draw(group) - if (returnRobj == FALSE) { + if (!returnRobj) { grDevices::pdf(outfile_Name, width = 23.4, height = 16.5) grid.draw(group) grDevices::dev.off() diff --git a/R/forest_plot2.R b/R/forest_plot2.R index 90200c00..645b5e7c 100644 --- a/R/forest_plot2.R +++ b/R/forest_plot2.R @@ -242,7 +242,7 @@ forest_plot_basic <- function( # OR or log(OR)? # If CI are symmetric then log(OR) # Use this to guess where to put the null line - null_line <- ifelse(all.equal(dat$effect - dat$lo_ci, dat$up_ci - dat$effect) == TRUE, 0, 1) + null_line <- if (isTRUE(all.equal(dat$effect - dat$lo_ci, dat$up_ci - dat$effect))) 0 else 1 # Change lab if (!is.null(xlim)) { @@ -435,7 +435,7 @@ forest_plot_names <- function(dat, section = NULL, bottom = TRUE) { # OR or log(OR)? # If CI are symmetric then log(OR) # Use this to guess where to put the null line - null_line <- ifelse(all.equal(dat$effect - dat$lo_ci, dat$up_ci - dat$effect) == TRUE, 0, 1) + null_line <- if (isTRUE(all.equal(dat$effect - dat$lo_ci, dat$up_ci - dat$effect))) 0 else 1 # up <- max(dat$up_ci, na.rm=TRUE) # lo <- min(dat$lo_ci, na.rm=TRUE) diff --git a/R/forest_plot_1-to-many.R b/R/forest_plot_1-to-many.R index c4e98568..18918059 100644 --- a/R/forest_plot_1-to-many.R +++ b/R/forest_plot_1-to-many.R @@ -213,9 +213,9 @@ sort_1_to_many <- function( ) Letters <- sort(c(paste0("A", Letters), paste0("B", Letters), paste0("C", Letters))) groups <- unique(mr_res[, group]) - mr_res$Index <- unlist(lapply(seq_along(unique(mr_res[, group])), FUN = function(x) { - rep(Letters[Letters == Letters[x]], length(which(mr_res[, group] == groups[x]))) - })) + # Assign each row the letter code for its group; match() makes this robust + # to groups whose rows are not contiguous in mr_res. + mr_res$Index <- Letters[match(mr_res[, group], groups)] mr_res <- mr_res[order(mr_res[, b], decreasing = TRUE), ] mr_res$Index2 <- Letters[seq_len(nrow(mr_res))] mr_res$Index3 <- paste(mr_res$Index, mr_res$Index2, sep = "") @@ -329,7 +329,7 @@ forest_plot_basic2 <- function( # OR or log(OR)? # If CI are symmetric then log(OR) # Use this to guess where to put the null line - null_line <- ifelse(all.equal(dat$effect - dat$lo_ci, dat$up_ci - dat$effect) == TRUE, 0, 1) + null_line <- if (isTRUE(all.equal(dat$effect - dat$lo_ci, dat$up_ci - dat$effect))) 0 else 1 # Change lab if (!is.null(xlim)) { @@ -565,7 +565,7 @@ forest_plot_names2 <- function( # OR or log(OR)? # If CI are symmetric then log(OR) # Use this to guess where to put the null line - null_line <- ifelse(all.equal(dat$effect - dat$lo_ci, dat$up_ci - dat$effect) == TRUE, 0, 1) + null_line <- if (isTRUE(all.equal(dat$effect - dat$lo_ci, dat$up_ci - dat$effect))) 0 else 1 # up <- max(dat$up_ci, na.rm=TRUE) # lo <- min(dat$lo_ci, na.rm=TRUE) @@ -683,7 +683,7 @@ forest_plot_addcol <- function( # OR or log(OR)? # If CI are symmetric then log(OR) # Use this to guess where to put the null line - null_line <- ifelse(all.equal(dat$effect - dat$lo_ci, dat$up_ci - dat$effect) == TRUE, 0, 1) + null_line <- if (isTRUE(all.equal(dat$effect - dat$lo_ci, dat$up_ci - dat$effect))) 0 else 1 lo <- 0 up <- 1 diff --git a/R/format_mr_results2.R b/R/format_mr_results2.R index 683e0457..1ca73613 100644 --- a/R/format_mr_results2.R +++ b/R/format_mr_results2.R @@ -306,7 +306,7 @@ power_prune <- function(dat, method = 1, dist.outcome = "binary") { if (is.null(ncase)) { ncase <- NA } - if (any(is.na(ncase))) { + if (anyNA(ncase)) { ncase <- dat1$samplesize.outcome if (dist.outcome == "binary") { warning(paste( @@ -314,18 +314,22 @@ power_prune <- function(dat, method = 1, dist.outcome = "binary") { )) } } - if (any(is.na(ncase))) { + if (anyNA(ncase)) { stop("sample size missing for at least 1 summary set") } - dat1 <- dat1[order(ncase, decreasing = TRUE), ] + # Permute dat1 and ncase identically; sort() would drop NA values and + # misalign ncase from dat1's rows. + ncase_ord <- order(ncase, decreasing = TRUE) + dat1 <- dat1[ncase_ord, ] # id.expout<-paste(split_exposure(dat)$exposure,split_outcome(dat)$outcome) - ncase <- ncase[order(ncase, decreasing = TRUE)] + ncase <- ncase[ncase_ord] # dat1$power.prune.ncase<-"drop" # dat1$power.prune.ncase[ncase==ncase[1]]<-"keep" dat1 <- dat1[ncase == ncase[1], ] nexp <- dat1$samplesize.exposure - dat1 <- dat1[order(nexp, decreasing = TRUE), ] - nexp <- nexp[order(nexp, decreasing = TRUE)] + nexp_ord <- order(nexp, decreasing = TRUE) + dat1 <- dat1[nexp_ord, ] + nexp <- nexp[nexp_ord] # dat1$power.prune.nexp<-"drop" # dat1$power.prune.nexp[nexp==nexp[1]]<-"keep" # dat1$power.prune<-"drop" @@ -365,14 +369,14 @@ power_prune <- function(dat, method = 1, dist.outcome = "binary") { z <- dat2$beta.exposure / dat2$se.exposure n <- dat2$samplesize.exposure b <- z / sqrt(2 * p * (1 - p) * (n + z^2)) - if (any(is.na(dat2$ncase.outcome))) { + if (anyNA(dat2$ncase.outcome)) { stop(paste("number of cases missing for summary set: ", id.subset.unique[j], sep = "")) } n.cas <- dat2$ncase.outcome n.con <- dat2$ncontrol.outcome var <- 1 # variance of risk factor assumed to be 1 r2 <- 2 * b^2 * p * (1 - p) / var - if (any(is.na(r2))) { + if (anyNA(r2)) { warning( "beta or allele frequency missing for some SNPs, which could affect accuracy of power pruning" ) @@ -386,7 +390,7 @@ power_prune <- function(dat, method = 1, dist.outcome = "binary") { iv.se <- 1 / sqrt(mean(dat2$samplesize.outcome, na.rm = TRUE) * r2sum) #standard error of the IV should be proportional to this } if (dist.outcome == "binary") { - if (any(is.na(n.cas)) || any(is.na(n.con))) { + if (anyNA(n.cas) || anyNA(n.con)) { warning( "dist.outcome set to binary but number of cases or controls is missing. Will try using total sample size instead but power pruning will be less accurate" ) diff --git a/R/ldsc.r b/R/ldsc.r index d6841bb1..a781cfa3 100644 --- a/R/ldsc.r +++ b/R/ldsc.r @@ -191,7 +191,7 @@ ldsc_rg <- function(id1, id2, ancestry = "infer", snpinfo = NULL, splitsize = 20 dplyr::inner_join(., d2, by = "rsid") h1 <- ldsc_h2_internal(d1$z1, d1$l2, d1$n1, d1$l2) - h2 <- ldsc_h2_internal(d2$z2, d2$l2, d2$n2, d1$l2) + h2 <- ldsc_h2_internal(d2$z2, d2$l2, d2$n2, d2$l2) dat <- dplyr::inner_join(d1, d2, by = "rsid") %>% dplyr::mutate( @@ -222,8 +222,9 @@ ldsc_rg <- function(id1, id2, ancestry = "infer", snpinfo = NULL, splitsize = 20 extract_split <- function(snplist, id, splitsize = 20000) { - nsplit <- round(length(snplist) / splitsize) - split(snplist, 1:nsplit) %>% + n <- length(snplist) + chunk_id <- rep(seq_len(ceiling(n / splitsize)), each = splitsize)[seq_len(n)] + split(snplist, chunk_id) %>% pbapply::pblapply(., function(x) { ieugwasr::associations(x, id, proxies = FALSE, x_api_source = x_api_source_header()) }) %>% diff --git a/R/leaveoneout.R b/R/leaveoneout.R index 706a048b..bc5f23f0 100644 --- a/R/leaveoneout.R +++ b/R/leaveoneout.R @@ -25,7 +25,7 @@ mr_leaveoneout <- function(dat, parameters = default_parameters(), method = mr_i exp_id <- combos$id.exposure[i] out_id <- combos$id.outcome[i] X <- dat_dt[id.exposure == exp_id & id.outcome == out_id] - x <- X[mr_keep == TRUE] + x <- X[X$mr_keep] nsnp <- nrow(x) if (nsnp == 0) { x <- X[1, ] diff --git a/R/mr.R b/R/mr.R index 7148d70c..f9c27f9b 100644 --- a/R/mr.R +++ b/R/mr.R @@ -58,7 +58,7 @@ mr <- function( exp_id <- combos$id.exposure[i] out_id <- combos$id.outcome[i] x1 <- dat_dt[id.exposure == exp_id & id.outcome == out_id] - x <- x1[mr_keep == TRUE] + x <- x1[x1$mr_keep] if (nrow(x) == 0) { message( @@ -818,8 +818,11 @@ mr_simple_median <- function(b_exp, b_out, se_exp, se_out, parameters = default_ #' @export #' @return MR estimate weighted_median <- function(b_iv, weights) { - betaIV.order <- b_iv[order(b_iv)] - weights.order <- weights[order(b_iv)] + # Order b_iv and weights by the same permutation so they stay aligned; sort() + # would drop NA b_iv values and misalign the two vectors. + ord <- order(b_iv) + betaIV.order <- b_iv[ord] + weights.order <- weights[ord] weights.sum <- cumsum(weights.order) - 0.5 * weights.order weights.sum <- weights.sum / sum(weights.order) below <- max(which(weights.sum < 0.5)) diff --git a/R/mr_mode.R b/R/mr_mode.R index 344841d0..6f6d2988 100644 --- a/R/mr_mode.R +++ b/R/mr_mode.R @@ -95,17 +95,20 @@ mr_mode <- function(dat, parameters = default_parameters(), mode_method = "all") phi = phi ) #Penalised mode, not assuming NOME - penalty <- stats::pchisq( - weights * (BetaIV.boot - beta.boot[i, (nphi + 1):(2 * nphi)])^2, - df = 1, - lower.tail = FALSE - ) - pen.weights <- weights * pmin(1, penalty * parameters$penk) - - beta.boot[i, (2 * nphi + 1):(3 * nphi)] <- beta( - BetaIV.in = BetaIV.boot, - seBetaIV.in = sqrt(1 / pen.weights), - phi = phi + #Penalise each SNP against the weighted-mode estimate for the same phi + #value, so the per-SNP penalty is not recycled across phi when nphi > 1. + beta.boot[i, (2 * nphi + 1):(3 * nphi)] <- vapply( + seq_len(nphi), + function(j) { + penalty <- stats::pchisq( + weights * (BetaIV.boot - beta.boot[i, nphi + j])^2, + df = 1, + lower.tail = FALSE + ) + pen.weights <- weights * pmin(1, penalty * parameters$penk) + beta(BetaIV.in = BetaIV.boot, seBetaIV.in = sqrt(1 / pen.weights), phi = phi[j]) + }, + numeric(1) ) #Simple mode, assuming NOME @@ -144,10 +147,21 @@ mr_mode <- function(dat, parameters = default_parameters(), mode_method = "all") #Point causal effect estimate using the weighted mode (not asusming NOME) beta_WeightedMode <- beta(BetaIV.in = BetaIV, seBetaIV.in = seBetaIV[, 1], phi = phi) weights <- 1 / seBetaIV[, 1]^2 - penalty <- stats::pchisq(weights * (BetaIV - beta_WeightedMode)^2, df = 1, lower.tail = FALSE) - pen.weights <- weights * pmin(1, penalty * parameters$penk) # penalized - - beta_PenalisedMode <- beta(BetaIV.in = BetaIV, seBetaIV.in = sqrt(1 / pen.weights), phi = phi) + # Penalise each SNP against the weighted-mode estimate for the same phi value, + # so the per-SNP penalty is not recycled across phi when length(phi) > 1. + beta_PenalisedMode <- vapply( + seq_along(phi), + function(j) { + penalty <- stats::pchisq( + weights * (BetaIV - beta_WeightedMode[j])^2, + df = 1, + lower.tail = FALSE + ) + pen.weights <- weights * pmin(1, penalty * parameters$penk) # penalized + beta(BetaIV.in = BetaIV, seBetaIV.in = sqrt(1 / pen.weights), phi = phi[j]) + }, + numeric(1) + ) #Point causal effect estimate using the weighted mode (asusming NOME) beta_WeightedMode_NOME <- beta(BetaIV.in = BetaIV, seBetaIV.in = seBetaIV[, 2], phi = phi) diff --git a/R/other_formats.R b/R/other_formats.R index ab48b67c..1c3a2122 100644 --- a/R/other_formats.R +++ b/R/other_formats.R @@ -104,7 +104,7 @@ harmonise_ld_dat <- function(x, ld) { return(NULL) } - if (any(!snpnames$keep)) { + if (!all(snpnames$keep)) { message( " - the following SNPs could not be aligned to the LD reference panel: \n- ", paste(subset(snpnames, !keep)$SNP, collapse = "\n - ") @@ -123,7 +123,7 @@ harmonise_ld_dat <- function(x, ld) { rownames(ld) <- snpnames$SNP colnames(ld) <- snpnames$SNP - if (any(!snpnames$keep)) { + if (!all(snpnames$keep)) { message("Removing ", sum(!snpnames$keep), " variants due to harmonisation issues") ld <- ld[snpnames$keep, snpnames$keep] x <- x[snpnames$keep, ] diff --git a/R/query.R b/R/query.R index 2f7471a2..36cc3f2d 100644 --- a/R/query.R +++ b/R/query.R @@ -115,18 +115,15 @@ extract_outcome_data_internal <- function( ) outcomes <- unique(outcomes) - if (proxies == FALSE) { + if (proxies %in% c(FALSE, 0)) { proxies <- 0 - } else if (proxies == TRUE) { + } else if (proxies %in% c(TRUE, 1)) { proxies <- 1 } else { stop("'proxies' argument should be TRUE or FALSE") } - if ( - (length(snps) < splitsize && length(outcomes) < splitsize) || - (length(outcomes) < splitsize && length(snps) < splitsize) - ) { + if (length(snps) < splitsize && length(outcomes) < splitsize) { d <- ieugwasr::associations( variants = snps, id = outcomes, @@ -318,14 +315,14 @@ format_d <- function(d) { mrcols <- c("beta.outcome", "se.outcome", "effect_allele.outcome") d$mr_keep.outcome <- complete.cases(d[, mrcols]) - if (any(!d$mr_keep.outcome)) { + if (!all(d$mr_keep.outcome)) { missinginfosnps <- paste(subset(d, !mr_keep.outcome)$SNP, collapse = " ") warning( "The following SNP(s) are missing required information for the MR tests and will be excluded: ", missinginfosnps ) } - if (all(!d$mr_keep.outcome)) { + if (!any(d$mr_keep.outcome)) { warning( "None of the provided SNPs can be used for MR analysis, they are missing required information." ) diff --git a/R/read_data.R b/R/read_data.R index fa64fe1e..628930f4 100644 --- a/R/read_data.R +++ b/R/read_data.R @@ -448,7 +448,7 @@ format_data <- function( dat$pval.outcome[index] <- min_pval dat$pval_origin.outcome <- "reported" - if (any(is.na(dat$pval.outcome))) { + if (anyNA(dat$pval.outcome)) { if ("beta.outcome" %in% names(dat) && "se.outcome" %in% names(dat)) { index <- is.na(dat$pval.outcome) dat$pval.outcome[index] <- 2 * @@ -551,7 +551,7 @@ format_data <- function( mrcols_present <- mrcols[mrcols %in% names(dat)] dat$mr_keep.outcome <- dat$mr_keep.outcome & complete.cases(dat[, mrcols_present]) - if (any(!dat$mr_keep.outcome)) { + if (!all(dat$mr_keep.outcome)) { missinginfosnps <- paste(subset(dat, !mr_keep.outcome)$SNP, collapse = " ") warning( "The following SNP(s) are missing required information for the MR tests and will be excluded: ", @@ -559,7 +559,7 @@ format_data <- function( ) } } - if (all(!dat$mr_keep.outcome)) { + if (!any(dat$mr_keep.outcome)) { warning( "None of the provided SNPs can be used for MR analysis, they are missing required information." ) @@ -748,7 +748,7 @@ random_string <- function(n = 1, len = 6) { create_ids <- function(x) { a <- as.factor(x) - levels(a) <- random_string(length(levels(a))) + levels(a) <- random_string(nlevels(a)) a <- as.character(a) return(a) } @@ -775,12 +775,12 @@ combine_data <- function(x) { stop("Datasets must be generated from format_data") } - check <- all(sapply(x, function(i) { + check <- all(vapply(x, function(i) { type %in% names(i) - })) + }, logical(1))) if (!check) { - stop("Not all datasets or of type '", type, "'") + stop("Not all datasets are of type '", type, "'") } id_col <- paste0("id.", type) diff --git a/R/rucker.R b/R/rucker.R index d964b6a6..b9970a47 100644 --- a/R/rucker.R +++ b/R/rucker.R @@ -40,9 +40,7 @@ PM <- function(y = y, s = s, Alpha = 0.1) { for (j in seq_len(L)) { tausq <- 0 Fstat <- 1 - TAUsq <- NULL while (Fstat > 0) { - TAUsq <- c(TAUsq, tausq) w <- 1 / (s^2 + tausq) sum.w <- sum(w) w2 <- w^2 diff --git a/R/scatterplot.R b/R/scatterplot.R index abc6dc88..f70a78d6 100644 --- a/R/scatterplot.R +++ b/R/scatterplot.R @@ -11,7 +11,7 @@ mr_scatter_plot <- function(mr_results, dat) { combos <- unique(dat[, c("id.exposure", "id.outcome")]) mrres <- lapply(seq_len(nrow(combos)), function(i) { d <- dat[dat$id.exposure == combos$id.exposure[i] & dat$id.outcome == combos$id.outcome[i], ] - if (nrow(d) < 2 | sum(d$mr_keep) == 0) { + if (nrow(d) < 2 || sum(d$mr_keep) == 0) { return(blank_plot("Insufficient number of SNPs")) } d <- subset(d, mr_keep) diff --git a/R/singlesnp.R b/R/singlesnp.R index 23ea9622..08a54707 100644 --- a/R/singlesnp.R +++ b/R/singlesnp.R @@ -41,7 +41,7 @@ mr_singlesnp <- function( exp_id <- combos$id.exposure[i] out_id <- combos$id.outcome[i] X <- dat_dt[id.exposure == exp_id & id.outcome == out_id] - x <- X[mr_keep == TRUE] + x <- X[X$mr_keep] nsnp <- nrow(x) if (nsnp == 0) { x <- X[1, ] diff --git a/R/steiger.R b/R/steiger.R index 429311f9..d2505066 100644 --- a/R/steiger.R +++ b/R/steiger.R @@ -249,7 +249,7 @@ directionality_test <- function(dat) { #' \item{sensitivity_plot}{Plot of parameter space of causal directions and measurement error} #' } mr_steiger2 <- function(r_exp, r_out, n_exp, n_out, r_xxo = 1, r_yyo = 1, ...) { - index <- any(is.na(r_exp)) | any(is.na(r_out)) | any(is.na(n_exp)) | any(is.na(n_out)) + index <- anyNA(r_exp) | anyNA(r_out) | anyNA(n_exp) | anyNA(n_out) n_exp <- n_exp[!index] n_out <- n_out[!index] diff --git a/tests/testthat/test_add_metadata.r b/tests/testthat/test_add_metadata.r index 4eb0e0ab..2acd4777 100644 --- a/tests/testthat/test_add_metadata.r +++ b/tests/testthat/test_add_metadata.r @@ -24,7 +24,7 @@ load(system.file("extdata", "test_add_metadata.RData", package = "TwoSampleMR")) test_that("exposure data 1", { d1 <- try(d1 %>% add_metadata()) - if (class(d1) == "try-error") { + if (inherits(d1, "try-error")) { skip("Server issues") } expect_true("units.exposure" %in% names(d1)) @@ -32,7 +32,7 @@ test_that("exposure data 1", { test_that("exposure data 2", { d2 <- try(d2 %>% add_metadata()) - if (class(d2) == "try-error") { + if (inherits(d2, "try-error")) { skip("Server issues") } expect_true("units.exposure" %in% names(d2)) @@ -40,7 +40,7 @@ test_that("exposure data 2", { test_that("outcome data 1", { d3 <- try(d3 %>% add_metadata()) - if (class(d3) == "try-error") { + if (inherits(d3, "try-error")) { skip("Server issues") } expect_true("units.outcome" %in% names(d3)) @@ -48,7 +48,7 @@ test_that("outcome data 1", { test_that("outcome data 2", { d4 <- try(d4 %>% add_metadata()) - if (class(d4) == "try-error") { + if (inherits(d4, "try-error")) { skip("Server issues") } expect_true("units.outcome" %in% names(d4)) @@ -56,7 +56,7 @@ test_that("outcome data 2", { test_that("dat 2", { d5 <- try(d5 %>% add_metadata()) - if (class(d5) == "try-error") { + if (inherits(d5, "try-error")) { skip("Server issues") } expect_true("units.outcome" %in% names(d5) & "units.exposure" %in% names(d5)) @@ -65,7 +65,7 @@ test_that("dat 2", { test_that("no id1", { d6$id.exposure <- "not a real id" d6 <- try(add_metadata(d6)) - if (class(d6) == "try-error") { + if (inherits(d6, "try-error")) { skip("Server issues") } expect_true(!"units.exposure" %in% names(d6)) @@ -74,7 +74,7 @@ test_that("no id1", { test_that("no id2", { d7$id.outcome <- "not a real id" d7 <- try(add_metadata(d7)) - if (class(d7) == "try-error") { + if (inherits(d7, "try-error")) { skip("Server issues") } expect_true(!"units.outcome" %in% names(d7)) @@ -82,7 +82,7 @@ test_that("no id2", { test_that("ukb-d", { d8 <- try(add_metadata(d8)) - if (class(d8) == "try-error") { + if (inherits(d8, "try-error")) { skip("Server issues") } expect_true("units.outcome" %in% names(d8)) @@ -90,18 +90,18 @@ test_that("ukb-d", { test_that("bbj-a-1", { d9 <- try(d9 %>% add_metadata()) - if (class(d9) == "try-error") { + if (inherits(d9, "try-error")) { skip("Server issues") } expect_true("samplesize.exposure" %in% names(d9)) - expect_true(all(!is.na(d9$samplesize.exposure))) + expect_true(!anyNA(d9$samplesize.exposure)) }) test_that("ieu-b-109", { d10 <- try(d10 %>% add_metadata()) - if (class(d10) == "try-error") { + if (inherits(d10, "try-error")) { skip("Server issues") } expect_true("samplesize.exposure" %in% names(d10)) - expect_true(all(!is.na(d10$samplesize.exposure))) + expect_true(!anyNA(d10$samplesize.exposure)) }) diff --git a/tests/testthat/test_instruments.R b/tests/testthat/test_instruments.R index 4264692e..9509a201 100644 --- a/tests/testthat/test_instruments.R +++ b/tests/testthat/test_instruments.R @@ -8,7 +8,7 @@ skip_on_cran() test_that("server and mrinstruments 1", { # no no exp_dat <- try(extract_instruments(outcomes = c("ieu-a-1032"))) - if (class(exp_dat) == "try-error") { + if (inherits(exp_dat, "try-error")) { skip("Server issues") } expect_true(length(unique(exp_dat$id)) == 0) @@ -18,7 +18,7 @@ test_that("server and mrinstruments 1", { test_that("server and mrinstruments 2", { # no yes exp_dat <- try(extract_instruments(outcomes = c("ebi-a-GCST004634"))) - if (class(exp_dat) == "try-error") { + if (inherits(exp_dat, "try-error")) { skip("Server issues") } expect_true(length(unique(exp_dat$id)) == 1) @@ -27,7 +27,7 @@ test_that("server and mrinstruments 2", { test_that("server and mrinstruments 3", { # yes no exp_dat <- try(extract_instruments(outcomes = c("ieu-a-2", "ieu-a-1032"))) - if (class(exp_dat) == "try-error") { + if (inherits(exp_dat, "try-error")) { skip("Server issues") } expect_true(length(unique(exp_dat$id)) == 1) @@ -36,7 +36,7 @@ test_that("server and mrinstruments 3", { test_that("server and mrinstruments 4", { # yes yes exp_dat <- try(extract_instruments(outcomes = c("ieu-a-2", "ebi-a-GCST004634"))) - if (class(exp_dat) == "try-error") { + if (inherits(exp_dat, "try-error")) { skip("Server issues") } expect_true(length(unique(exp_dat$id)) == 2) @@ -44,7 +44,7 @@ test_that("server and mrinstruments 4", { test_that("server and mrinstruments 5", { exp_dat <- try(extract_instruments(outcomes = c("ieu-a-1032", "ebi-a-GCST004634"))) - if (class(exp_dat) == "try-error") { + if (inherits(exp_dat, "try-error")) { skip("Server issues") } expect_true(length(unique(exp_dat$id)) == 1) @@ -52,7 +52,7 @@ test_that("server and mrinstruments 5", { test_that("server and mrinstruments 6", { exp_dat <- try(extract_instruments(outcomes = c(2, 100, "ieu-a-1032", 104, 72, 999))) - if (class(exp_dat) == "try-error") { + if (inherits(exp_dat, "try-error")) { skip("Server issues") } expect_true(length(unique(exp_dat$id)) == 5) @@ -62,7 +62,7 @@ test_that("server and mrinstruments 7", { exp_dat <- try(extract_instruments( outcomes = c(2, 100, "ieu-a-1032", 104, 72, 999, "ebi-a-GCST004634") )) - if (class(exp_dat) == "try-error") { + if (inherits(exp_dat, "try-error")) { skip("Server issues") } expect_true(length(unique(exp_dat$id)) == 6) @@ -72,7 +72,7 @@ load(system.file("extdata", "test_commondata.RData", package = "TwoSampleMR")) test_that("read data", { exp_dat <- try(extract_instruments("ieu-a-2")) - if (class(exp_dat) == "try-error") { + if (inherits(exp_dat, "try-error")) { skip("Server issues") } names(exp_dat) <- gsub(".exposure", "", names(exp_dat)) diff --git a/tests/testthat/test_ld.R b/tests/testthat/test_ld.R index fa44a46e..4ed44bf3 100644 --- a/tests/testthat/test_ld.R +++ b/tests/testthat/test_ld.R @@ -7,11 +7,11 @@ skip_on_cran() # extract some data a <- try(extract_instruments("ieu-a-2", clump = 0)) -if (class(a) == "try-error") { +if (inherits(a, "try-error")) { skip("Server issues") } out <- try(clump_data(a)) -if (class(out) == "try-error") { +if (inherits(out, "try-error")) { skip("Server issues") } @@ -38,7 +38,7 @@ test_that("matrix", { test_that("clump multiple", { a <- try(extract_instruments(c("ieu-a-2", "ieu-a-1001"), clump = 0)) - if (class(a) == "try-error") { + if (inherits(a, "try-error")) { skip("Server issues") } out <- clump_data(a) diff --git a/tests/testthat/test_plots.R b/tests/testthat/test_plots.R index 00a1fd34..3255167e 100644 --- a/tests/testthat/test_plots.R +++ b/tests/testthat/test_plots.R @@ -5,7 +5,7 @@ load(system.file("extdata", "test_commondata.RData", package = "TwoSampleMR")) test_that("MR scatter plot for mr_ivw", { # dat <- make_dat("ieu-a-2", "ieu-a-7") m <- mr(dat, method_list = "mr_ivw") - expect_no_error(p <- mr_scatter_plot(m, dat)) + p <- expect_no_error(mr_scatter_plot(m, dat)) expect_true(is.list(p)) expect_true(length(p) == 1L) }) @@ -13,21 +13,21 @@ test_that("MR scatter plot for mr_ivw", { test_that("Scatter plot for default set of estimates", { # dat <- make_dat("ieu-a-2", "ieu-a-7") m2 <- mr(dat) - expect_no_error(p2 <- mr_scatter_plot(m2, dat)) + p2 <- expect_no_error(mr_scatter_plot(m2, dat)) expect_true(is.list(p2)) expect_true(length(p2) == 1L) }) test_that("Scatter plot for mr_grip", { m3 <- mr(dat, method_list = "mr_grip") - expect_no_error(p3 <- mr_scatter_plot(m3, dat)) + p3 <- expect_no_error(mr_scatter_plot(m3, dat)) expect_true(is.list(p3)) expect_true(length(p3) == 1L) }) test_that("A second scatter plot for mr_grip", { m4 <- mr(dat2, method_list = "mr_grip") - expect_no_error(p4 <- mr_scatter_plot(m4, dat2)) + p4 <- expect_no_error(mr_scatter_plot(m4, dat2)) expect_true(is.list(p4)) expect_true(length(p4) == 4L) }) @@ -36,17 +36,17 @@ res_single <- mr_singlesnp(dat) res_loo <- mr_leaveoneout(dat) test_that("Forest plot", { - expect_no_error(p5 <- mr_forest_plot(res_single)) + p5 <- expect_no_error(mr_forest_plot(res_single)) expect_true(length(p5) == 1L) }) test_that("Leave one out plot", { - expect_no_error(p6 <- mr_leaveoneout_plot(res_loo)) + p6 <- expect_no_error(mr_leaveoneout_plot(res_loo)) expect_true(length(p6) == 1L) }) test_that("Funnel plot", { - expect_no_error(p7 <- mr_funnel_plot(res_single)) + p7 <- expect_no_error(mr_funnel_plot(res_single)) expect_true(length(p7) == 1L) }) @@ -78,7 +78,7 @@ test_that("Forest plot 1 to many", { test_that("Forest plot 1 to many test 2", { res$pval <- formatC(res$pval, format = "e", digits = 2) expect_warning(utils::capture.output( - p9 <- forest_plot_1_to_many( + forest_plot_1_to_many( res, b = "b", se = "se", diff --git a/tests/testthat/test_rsq.r b/tests/testthat/test_rsq.r index bef89585..511e04b6 100644 --- a/tests/testthat/test_rsq.r +++ b/tests/testthat/test_rsq.r @@ -130,7 +130,7 @@ test_that("bbj-a-1", { if (inherits(d, "try-error")) { skip("Server issues") } - expect_true(all(!is.na(d$rsq.exposure))) + expect_true(!anyNA(d$rsq.exposure)) }) test_that("bsen vs pn", { diff --git a/vignettes/perform_mr.Rmd b/vignettes/perform_mr.Rmd index 637b0485..a7b94fdd 100644 --- a/vignettes/perform_mr.Rmd +++ b/vignettes/perform_mr.Rmd @@ -868,7 +868,7 @@ Alternatively, convert to the `MRInput` format but also obtaining the LD matrix ```{r eval=evalinr44} dat2 <- try(dat_to_MRInput(dat, get_correlation = TRUE)) -if (class(dat2) != "try-error") MendelianRandomization::mr_ivw(dat2[[1]], correl = TRUE) +if (!inherits(dat2, "try-error")) MendelianRandomization::mr_ivw(dat2[[1]], correl = TRUE) ``` ## MR-MoE: Using a mixture of experts machine learning approach