added mr_sign test

R/mr.R

@@ -187,13 +187,21 @@ mr_method_list <- function()
 			use_by_default=FALSE,
 			heterogeneity_test=FALSE
 		),
-                list(
+		list(
 			obj="mr_raps",
 			name="Robust adjusted profile score (RAPS)",
 			PubmedID="",
 			Description="",
 			use_by_default=FALSE,
 			heterogeneity_test=FALSE
+		),
+		list(
+			obj="mr_sign",
+			name="Sign concordance test",
+			PubmedID="",
+			Description="Tests for concordance of signs between exposure and outcome",
+			use_by_default=FALSE,
+			heterogeneity_test=FALSE
 		)
 	)
 	a <- lapply(a, as.data.frame)
@@ -885,3 +893,41 @@ mr_raps <- function(b_exp, b_out, se_exp, se_out, parameters = default_parameter
          nsnp = length(b_exp))
 
 }
+
+#' MR sign test
+#' 
+#' Tests how often the SNP-exposure and SNP-outcome signs are concordant
+#' This is to avoid the problem of averaging over all SNPs, which can suffer bias due to outliers with strong effects; and to avoid excluding SNPs which is implicit in median and mode based estimators
+#' The effect estimate here is not to be interpreted as the effect size - it is the proportion of SNP-exposure and SNP-outcome effects that have concordant signs. 
+#' e.g. +1 means all have the same sign, -1 means all have opposite signs, and 0 means that there is an equal number of concordant and discordant signs.
+#' Restricted to only work if there are 6 or more valud SNPs
+#' 
+#' @param b_exp Vector of genetic effects on exposure
+#' @param b_out Vector of genetic effects on outcome
+#' @param se_exp Not required
+#' @param se_out Not required
+#' @param parameters Not required
+#'
+#' @export
+#' @return List with the following elements:
+#'         b: Concordance (see description)
+#'         se: NA
+#'         pval: p-value
+#'         nsnp: Number of SNPs (excludes NAs and effect estimates that are 0)
+mr_sign <- function(b_exp, b_out, se_exp=NULL, se_out=NULL, parameters=NULL)
+{
+	b_exp[b_exp == 0] <- NA
+	b_out[b_out == 0] <- NA
+	if(sum(!is.na(b_exp) & !is.na(b_out)) < 6)
+	{
+		return(list(b=NA, se=NA, pval=NA, nsnp=NA))
+	}
+	x <- sum(sign(b_exp) == sign(b_out), na.rm=TRUE)
+	n <- sum(!is.na(b_exp) & !is.na(b_out))
+
+	out <- binom.test(x=x, n=n, p=0.5)
+	b <- (out$estimate - 0.5) * 2
+	names(b) <- NULL
+	pval <- out$p.value
+	return(list(b=b, se=NA, pval=pval, nsnp=n))
+}

tests/api_alive.r

@@ -9,3 +9,110 @@ d <- extract_outcome_data(a$SNP, 7)
 e <- available_outcomes()
 message(nrow(e))
 
+
+
+
+
+
+load("~/repo/mr-eve/results/03/exposure_dat.rdata")
+snplist <- unique(exposure_dat$SNP)
+
+t1 <- Sys.time()
+options(mrbaseapi="http://ieu-db-interface.epi.bris.ac.uk:8080/")
+out1 <- extract_outcome_data(snplist, 1001, proxies=FALSE)
+Sys.time()-t1
+
+
+t1 <- Sys.time()
+toggle_api("release")
+out2 <- extract_outcome_data(snplist, 1001, proxies=FALSE)
+Sys.time()-t1
+
+
+t1 <- Sys.time()
+options(mrbaseapi="http://ieu-db-interface.epi.bris.ac.uk:8080/")
+pout1 <- extract_outcome_data(snplist[1:5000], 1001, proxies=TRUE)
+Sys.time()-t1
+
+
+t1 <- Sys.time()
+toggle_api("release")
+pout2 <- extract_outcome_data(snplist[1:5000], 1001, proxies=TRUE)
+Sys.time()-t1
+
+
+
+library(dplyr)
+
+t1 <- Sys.time()
+options(mrbaseapi="http://ieu-db-interface.epi.bris.ac.uk:8080/")
+pout1 <- lapply(split(snplist[1:5000], 1:5), function(x) extract_outcome_data(x, 1001, proxies=TRUE)) %>% bind_rows()
+Sys.time()-t1
+# 48sec
+
+t1 <- Sys.time()
+options(mrbaseapi="http://ieu-db-interface.epi.bris.ac.uk:8080/")
+pout2 <- extract_outcome_data(snplist[1:5000], 1001, proxies=TRUE)
+Sys.time()-t1
+# 39sec
+
+t1 <- Sys.time()
+toggle_api("release")
+pout2 <- lapply(split(snplist[1:5000], 1:5), function(x) extract_outcome_data(x, 1001, proxies=TRUE)) %>% bind_rows()
+Sys.time()-t1
+# 53sec
+
+t1 <- Sys.time()
+toggle_api("release")
+pout2 <- extract_outcome_data(snplist[1:5000], 1001, proxies=TRUE)
+Sys.time()-t1
+# 40sec
+
+####
+
+t1 <- Sys.time()
+options(mrbaseapi="http://ieu-db-interface.epi.bris.ac.uk:8080/")
+pout1 <- lapply(split(snplist, 1:4), function(x) extract_outcome_data(x, 1001, proxies=TRUE)) %>% bind_rows()
+Sys.time()-t1
+# 3,7m
+
+t1 <- Sys.time()
+options(mrbaseapi="http://ieu-db-interface.epi.bris.ac.uk:8080/")
+pout2 <- extract_outcome_data(snplist, 1001, proxies=TRUE, proxy_splitsize=50)
+Sys.time()-t1
+# 2.6m - 500
+
+
+t1 <- Sys.time()
+options(mrbaseapi="http://ieu-db-interface.epi.bris.ac.uk:8080/")
+pout2 <- extract_outcome_data(snplist, 1001, proxies=TRUE, proxy_splitsize=100)
+Sys.time()-t1
+# 2.3m - 500
+
+
+
+t1 <- Sys.time()
+options(mrbaseapi="http://ieu-db-interface.epi.bris.ac.uk:8080/")
+pout2 <- extract_outcome_data(snplist, 1001, proxies=TRUE)
+Sys.time()-t1
+
+
+t1 <- Sys.time()
+toggle_api("release")
+pout2 <- extract_outcome_data(snplist, 1001, proxies=TRUE)
+Sys.time()-t1
+
+
+
+t1 <- Sys.time()
+toggle_api("release")
+pout2 <- lapply(split(snplist[1:5000], 1:5), function(x) extract_outcome_data(x, 1001, proxies=TRUE)) %>% bind_rows()
+Sys.time()-t1
+# 53sec
+
+t1 <- Sys.time()
+toggle_api("release")
+pout2 <- extract_outcome_data(snplist[1:5000], 1001, proxies=TRUE)
+Sys.time()-t1
+# 40sec
+

tests/test_mr_sign.R

@@ -0,0 +1,32 @@
+
+# MR sign concordance - how often are the SNP-exposure and SNP-outcome effects in the same direction?
+
+# Use binomial test assuming that on average should be 50% under the null hypothesis
+
+# What is the minumum number of SNPs required to achieve p-value of 0.05?
+
+param <- expand.grid(n = 1:100, x=0:100)
+param <- subset(param, x <= n)
+for(i in 1:nrow(param))
+{
+	param$pval[i] <- binom.test(x=param$x[i], n=param$n[i], p=0.5)$p.value
+}
+min(param$pval)
+
+library(ggplot2)
+ggplot(param, aes(x=x, y=n)) +
+geom_point(aes(colour=pval < 0.05))
+
+library(dplyr)
+group_by(param, n) %>% summarise(minp = min(pval))
+
+# Looks like 6 is the smallest number of SNPs that can achieve p < 0.05
+
+a <- extract_instruments(2)
+b <- extract_outcome_data(a$SNP,7)
+dat <- harmonise_data(a, b)
+mr(dat, method_list=c("mr_ivw", "mr_sign"))
+
+
+# Could consider doing a parametric bootstrap type analysis
+