diff --git a/DESCRIPTION b/DESCRIPTION
index 59f5c661..165043bd 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -47,7 +47,8 @@ Imports:
     psych,
     RadialMR,
     reshape2,
-    rmarkdown
+    rmarkdown,
+    rlang
 Suggests: 
     Cairo,
     car,
diff --git a/NAMESPACE b/NAMESPACE
index df469578..0b52b4e5 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -109,8 +109,12 @@ export(split_outcome)
 export(standardise_units)
 export(steiger_filtering)
 export(steiger_sensitivity)
+export(steiger_sensitivity_conf)
 export(subset_on_method)
 export(trim)
 export(weighted_median)
 export(weighted_median_bootstrap)
 importFrom(magrittr,"%>%")
+importFrom(rlang,.data)
+importFrom(stats,dbeta)
+importFrom(stats,quantile)
diff --git a/R/leaveoneout.R b/R/leaveoneout.R
index a598d6ca..080ca467 100644
--- a/R/leaveoneout.R
+++ b/R/leaveoneout.R
@@ -112,11 +112,11 @@ mr_leaveoneout_plot <- function(leaveoneout_results)
 		ggplot2::ggplot(d, ggplot2::aes(y=SNP, x=b)) +
 		ggplot2::geom_vline(xintercept=0, linetype="dotted") +
 		# ggplot2::geom_errorbarh(ggplot2::aes(xmin=pmax(lo, min(d$b, na.rm=T)), xmax=pmin(up, max(d$b, na.rm=T)), size=as.factor(tot), colour=as.factor(tot)), height=0) +
-		ggplot2::geom_errorbarh(ggplot2::aes(xmin=lo, xmax=up, size=as.factor(tot), colour=as.factor(tot)), height=0) +
+		ggplot2::geom_errorbarh(ggplot2::aes(xmin=lo, xmax=up, linewidth=as.factor(tot), colour=as.factor(tot)), height=0) +
 		ggplot2::geom_point(ggplot2::aes(colour=as.factor(tot))) +
 		ggplot2::geom_hline(ggplot2::aes(yintercept = which(levels(SNP) %in% "")), colour="grey") +
 		ggplot2::scale_colour_manual(values=c("black", "red")) +
-		ggplot2::scale_size_manual(values=c(0.3, 1)) +
+		ggplot2::scale_linewidth_manual(values=c(0.3, 1)) +
 		# xlim(c(min(c(0, d$b), na.rm=T), max(c(0, d$b), na.rm=T))) +
 		ggplot2::theme(
 			legend.position="none", 
diff --git a/R/other_formats.R b/R/other_formats.R
index adb0fda3..37fedc96 100644
--- a/R/other_formats.R
+++ b/R/other_formats.R
@@ -172,7 +172,7 @@ dat_to_RadialMR <- function(dat)
 		message("Converting:")
 		message(" - exposure: ", x$exposure[1])
 		message(" - outcome: ", x$outcome[1])
-		d <- subset(x, mr_keep=TRUE)
+		d <- subset(x, mr_keep)
 		d <- RadialMR::format_radial(d$beta.exposure, d$beta.outcome, d$se.exposure, d$se.outcome, RSID=d$SNP)
 		return(d)
 	})
diff --git a/R/singlesnp.R b/R/singlesnp.R
index 301d70a5..942ff53c 100644
--- a/R/singlesnp.R
+++ b/R/singlesnp.R
@@ -114,11 +114,11 @@ mr_forest_plot <- function(singlesnp_results, exponentiate=FALSE)
 		ggplot2::ggplot(d, ggplot2::aes(y=SNP, x=b)) +
 		ggplot2::geom_vline(xintercept=xint, linetype="dotted") +
 		# ggplot2::geom_errorbarh(ggplot2::aes(xmin=pmax(lo, min(d$b, na.rm=T)), xmax=pmin(up, max(d$b, na.rm=T)), size=as.factor(tot), colour=as.factor(tot)), height=0) +
-		ggplot2::geom_errorbarh(ggplot2::aes(xmin=lo, xmax=up, size=as.factor(tot), colour=as.factor(tot)), height=0) +
+		ggplot2::geom_errorbarh(ggplot2::aes(xmin=lo, xmax=up, linewidth=as.factor(tot), colour=as.factor(tot)), height=0) +
 		ggplot2::geom_point(ggplot2::aes(colour=as.factor(tot))) +
 		ggplot2::geom_hline(ggplot2::aes(yintercept = which(levels(SNP) %in% "")), colour="grey") +
 		ggplot2::scale_colour_manual(values=c("black", "red")) +
-		ggplot2::scale_size_manual(values=c(0.3, 1)) +
+		ggplot2::scale_linewidth_manual(values=c(0.3, 1)) +
 		# xlim(c(min(c(0, d$b), na.rm=T), max(c(0, d$b), na.rm=T))) +
 		ggplot2::theme(
 			legend.position="none", 
diff --git a/R/steiger.R b/R/steiger.R
index 663761e9..849151d5 100644
--- a/R/steiger.R
+++ b/R/steiger.R
@@ -64,6 +64,125 @@ steiger_sensitivity <- function(rgx_o, rgy_o, ...)
 	))
 }
 
+ss_conf_calcs <- function(bxy, bgx, bux, buy, vg, vu, vex, vey) {
+    args <- as.list(environment())
+    bxyo <- ((bgx^2*bxy*vg + bux^2*bxy*vu + bux*buy*vu + bxy*vex) / (vg*bgx^2 + vu*bux^2 + vex))
+    vx <- bgx^2 * vg + bux^2 * vu + vex
+    vy <- bxy^2*bgx^2*vg + (bxy*bux+buy)^2*vu + bxy^2*(vex) + vey
+    conf <- bux * vu * buy / (vg * bgx^2 + vu * bux^2 + vex)
+    rsqxyo <- bxyo^2 * vx / vy
+    rsqxyos <- rsqxyo * sign(bxyo)
+    rsqxy <- bxy^2 * vx / vy
+    rsqxys <- rsqxy * sign(bxy)
+    rsqgx <- bgx^2*vg / (bgx^2 * vg + bux^2 * vu + vex)
+    rsqgy <- bgx^2*bxy^2*vg / (bxy^2*bgx^2*vg + (bxy*bux+buy)^2*vu + bxy^2*(vex) + vey)
+    rsqux <- bux^2*vu / (bgx^2 * vg + bux^2 * vu + vex)
+    rsquy <- (buy + bux * bxy)^2 * vu / (bxy^2*bgx^2*vg + (bxy*bux+buy)^2*vu + bxy^2*(vex) + vey)
+    rsquxs <- rsqux * sign(bux)
+    rsquys <- rsquy * sign(buy)
+    return(c(args, list(
+        vx=vx,
+        vy=vy,
+        bxyo=bxyo, 
+        conf=conf, 
+        rsqgx=rsqgx, 
+        rsqgy=rsqgy, 
+        rsqux=rsqux, 
+        rsquy=rsquy, 
+        rsqxy=rsqxy, 
+        rsqxyo=rsqxyo, 
+        rsqxyos=rsqxyos, 
+        rsquxs=rsquxs, 
+        rsquys=rsquys, 
+        rsqxys=rsqxys
+    )))
+}
+
+ss_conf_1d <- function(bxy=0.1, bxyo=0.2, bgx=0.5, vx=1, vy=1, vu=1, vg=0.5, simsize=100) {
+    # vx <- bgx^2 * vg + p$bux_vec^2 * vu + vex
+    bux_lim <- sqrt((vx - bgx^2 * vg)/vu)
+    bux_vec <- seq(-bux_lim, bux_lim, length.out=simsize)
+    # Allow causal effect to vary by +/- 200%
+    vex <- vx - bgx^2 * vg - bux_vec^2 * vu
+    conf <- bxyo - bxy
+    buy_vec <- conf * (bgx^2*vg + bux_vec^2*vu + vex) / (bux_vec * vu)
+    vey <- vy - (bxy^2*bgx^2*vg + (bxy*bux_vec+buy_vec)^2*vu + bxy^2*vex)
+    # vy <- bxy^2*bgx^2*vg + (bxy*bux_vec+buy_vec)^2*vu + bxy^2*vex + vey
+    bux_vec * vu * buy_vec / (vg * bgx^2 + vu * bux_vec^2 + vex)
+    res <- ss_conf_calcs(bxy, bgx, bux_vec, buy_vec, vg, vu, vex, vey) %>% 
+      dplyr::as_tibble() %>%
+      dplyr::mutate(bxy=bxy, bgx=bgx, bux=bux_vec, buy=buy_vec, vg=vg, vu=vu, vex=vex, vey=vey)
+    return(res)
+}
+
+#' Sensitivity analysis for unmeasured confounding on Steiger inferred causal direction
+#' 
+#' @description 
+#' This function takes known parameters from an MR analysis and determines the range of unmeasured 
+#' confounding that would be required to agree with the Steiger test inference, and the range of 
+#' unmeasured confounding values that would be required to disagree with the Steiger test inference
+#' 
+#' @param bxy MR estimate of x -> y
+#' @param bxyo Observational estimate of x -> y, if not available can re-run with different values to test sensitivity
+#' @param bgx SNP-exposure association
+#' @param vx Variance of X
+#' @param vy Variance of Y
+#' @param vg Variance of SNP (approximately 2*p*(1-p) where p is the allele frequency)
+#' @param vu Arbitrary variance of unmeasured confounder, default = 1
+#' @param simsize Density of search grid, default=10000
+#' @param beta_a Weighting of confounder values, which are beta distributed. Specify 'a' parameter of beta distribution, default=1 implying flat prior
+#' @param beta_b Weighting of confounder values, which are beta distributed. Specify 'b' parameter of beta distribution, default=1 implying flat prior
+#' @param plot Whether to generate plot. Default=TRUE
+#' 
+#' @importFrom stats dbeta quantile
+#' @importFrom rlang .data
+#' @export
+#' @return List of results
+#' \describe{
+#' \item{o}{data frame of possible confounding values and agreement with inferred steiger result}
+#' \item{prop}{(weighted) fraction of confounding space that agrees with inferred steiger result}
+#' \item{pl}{plot}
+#' }
+steiger_sensitivity_conf <- function(bxy, bxyo, bgx, vx, vy, vg, vu = 1, simsize=10000, beta_a=1, beta_b=1, plot=TRUE)
+{
+    o <- dplyr::bind_rows(
+        ss_conf_1d(bxy=bxy, bxyo=bxyo, bgx=bgx, vx=vx, vy=vy, vu=vu, vg=vg, simsize=simsize) %>%
+            dplyr::mutate(direction="inferred"),
+        ss_conf_1d(bxy=1/bxy, bxyo=bxyo * vx / vy, bgx=bgx * bxy, vx=vy, vy=vx, vu=vu, vg=vg, simsize=simsize) %>%
+            dplyr::mutate(direction="reverse")
+    ) %>%
+    dplyr::filter(
+        .data$vex >= 0 & 
+        .data$vey >= 0 &
+        .data$rsquy >= 0 & .data$rsquy <= 1 &
+        .data$rsqux >= 0 & .data$rsqux <= 1 &
+        .data$rsqgx >= 0 & .data$rsqgx <= 1 &
+        .data$rsqgy >= 0 & .data$rsqgy <= 1
+    ) %>%
+        dplyr::group_by(.data$direction) %>%
+        dplyr::do({
+            x <- .
+            x1 <- x$rsqux[-1]
+            x2 <- x$rsqux[-length(x$rsqux)]
+            y1 <- x$rsquy[-1]
+            y2 <- x$rsquy[-length(x$rsquy)]
+            d <- sqrt((x1-x2)^2 + (y1-y2)^2)
+            d[d > quantile(d, na.rm=T, probs=0.99)*4] <- NA
+            x$d <- c(NA, d)
+            x$weight <- dbeta(x$rsqux, shape1=beta_a, shape2=beta_b) * dbeta(x$rsquy, shape1=beta_a, shape2=beta_b)
+            x
+        })
+    
+    w <- o$d * o$weight
+    w1 <- w[o$direction=="inferred"]
+    prop <- sum(w1, na.rm=T) / sum(w, na.rm=T)
+    ret <- list(result=o, prop=prop)
+    if(plot) {
+        ret$pl <- ggplot2::ggplot(o, ggplot2::aes(x=.data$rsquxs, y=.data$rsquys)) +
+        ggplot2::geom_point(ggplot2::aes(colour=.data$direction, size=.data$weight))
+    }
+    return(ret)
+}
 
 #' MR Steiger test of directionality
 #'
diff --git a/man/steiger_sensitivity_conf.Rd b/man/steiger_sensitivity_conf.Rd
new file mode 100644
index 00000000..c94bd2a9
--- /dev/null
+++ b/man/steiger_sensitivity_conf.Rd
@@ -0,0 +1,56 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/steiger.R
+\name{steiger_sensitivity_conf}
+\alias{steiger_sensitivity_conf}
+\title{Sensitivity analysis for unmeasured confounding on Steiger inferred causal direction}
+\usage{
+steiger_sensitivity_conf(
+  bxy,
+  bxyo,
+  bgx,
+  vx,
+  vy,
+  vg,
+  vu = 1,
+  simsize = 10000,
+  beta_a = 1,
+  beta_b = 1,
+  plot = TRUE
+)
+}
+\arguments{
+\item{bxy}{MR estimate of x -> y}
+
+\item{bxyo}{Observational estimate of x -> y, if not available can re-run with different values to test sensitivity}
+
+\item{bgx}{SNP-exposure association}
+
+\item{vx}{Variance of X}
+
+\item{vy}{Variance of Y}
+
+\item{vg}{Variance of SNP (approximately 2\emph{p}(1-p) where p is the allele frequency)}
+
+\item{vu}{Arbitrary variance of unmeasured confounder, default = 1}
+
+\item{simsize}{Density of search grid, default=10000}
+
+\item{beta_a}{Weighting of confounder values, which are beta distributed. Specify 'a' parameter of beta distribution, default=1 implying flat prior}
+
+\item{beta_b}{Weighting of confounder values, which are beta distributed. Specify 'b' parameter of beta distribution, default=1 implying flat prior}
+
+\item{plot}{Whether to generate plot. Default=TRUE}
+}
+\value{
+List of results
+\describe{
+\item{o}{data frame of possible confounding values and agreement with inferred steiger result}
+\item{prop}{(weighted) fraction of confounding space that agrees with inferred steiger result}
+\item{pl}{plot}
+}
+}
+\description{
+This function takes known parameters from an MR analysis and determines the range of unmeasured
+confounding that would be required to agree with the Steiger test inference, and the range of
+unmeasured confounding values that would be required to disagree with the Steiger test inference
+}
diff --git a/tests/testthat/test_steiger.R b/tests/testthat/test_steiger.R
index 5c995df1..ccb81190 100644
--- a/tests/testthat/test_steiger.R
+++ b/tests/testthat/test_steiger.R
@@ -20,3 +20,9 @@ test_that("steiger filtering", {
 	w <- steiger_filtering(w)
 	expect_true("steiger_pval" %in% names(w))
 })
+
+
+test_that("steiger sensitivity uconf", {
+	r <- steiger_sensitivity_conf(bxy=0.1, bxyo=0.2, bgx=0.5, vx=1, vy=1, vg=0.5, beta_a=1, beta_b=10)
+	expect_equal(r$prop, 1)
+})