|
1 | | -#' Create scatter plot with lines showing the causal estimate for different MR tests |
| 1 | +#' Create scatter plot with fitted lines showing the causal effect estimate for different MR estimators |
2 | 2 | #' |
3 | | -#' Requires dev version of ggplot2 |
| 3 | +#' Create scatter plot with fitted lines showing the causal effect estimate for different MR estimators. |
4 | 4 | #' |
5 | 5 | #' @param mr_results Output from [mr()]. |
6 | 6 | #' @param dat Output from [harmonise_data()]. |
7 | 7 | #' @export |
8 | 8 | #' @return List of plots |
9 | | -mr_scatter_plot <- function(mr_results, dat) |
10 | | -{ |
11 | | - # dat <- subset(dat, paste(id.outcome, id.exposure) %in% paste(mr_results$id.outcome, mr_results$id.exposure)) |
12 | | - mrres <- plyr::dlply(dat, c("id.exposure", "id.outcome"), function(d) |
13 | | - { |
14 | | - d <- plyr::mutate(d) |
15 | | - if(nrow(d) < 2 | sum(d$mr_keep) == 0) |
16 | | - { |
17 | | - return(blank_plot("Insufficient number of SNPs")) |
18 | | - } |
19 | | - d <- subset(d, mr_keep) |
20 | | - index <- d$beta.exposure < 0 |
21 | | - d$beta.exposure[index] <- d$beta.exposure[index] * -1 |
22 | | - d$beta.outcome[index] <- d$beta.outcome[index] * -1 |
23 | | - mrres <- subset(mr_results, id.exposure == d$id.exposure[1] & id.outcome == d$id.outcome[1]) |
24 | | - mrres$a <- 0 |
25 | | - if("MR Egger" %in% mrres$method) |
26 | | - { |
27 | | - temp <- mr_egger_regression(d$beta.exposure, d$beta.outcome, d$se.exposure, d$se.outcome, default_parameters()) |
28 | | - mrres$a[mrres$method == "MR Egger"] <- temp$b_i |
29 | | - } |
| 9 | +mr_scatter_plot <- function(mr_results, dat) { |
| 10 | + # dat <- subset(dat, paste(id.outcome, id.exposure) %in% paste(mr_results$id.outcome, mr_results$id.exposure)) |
| 11 | + mrres <- plyr::dlply(dat, c("id.exposure", "id.outcome"), function(d) { |
| 12 | + d <- plyr::mutate(d) |
| 13 | + if (nrow(d) < 2 | sum(d$mr_keep) == 0) { |
| 14 | + return(blank_plot("Insufficient number of SNPs")) |
| 15 | + } |
| 16 | + d <- subset(d, mr_keep) |
| 17 | + index <- d$beta.exposure < 0 |
| 18 | + d$beta.exposure[index] <- d$beta.exposure[index] * -1 |
| 19 | + d$beta.outcome[index] <- d$beta.outcome[index] * -1 |
| 20 | + mrres <- subset( |
| 21 | + mr_results, |
| 22 | + id.exposure == d$id.exposure[1] & id.outcome == d$id.outcome[1] |
| 23 | + ) |
| 24 | + mrres$a <- 0 |
| 25 | + if ("MR Egger" %in% mrres$method) { |
| 26 | + temp <- mr_egger_regression( |
| 27 | + d$beta.exposure, |
| 28 | + d$beta.outcome, |
| 29 | + d$se.exposure, |
| 30 | + d$se.outcome, |
| 31 | + default_parameters() |
| 32 | + ) |
| 33 | + mrres$a[mrres$method == "MR Egger"] <- temp$b_i |
| 34 | + } |
30 | 35 |
|
31 | | - if("MR Egger (bootstrap)" %in% mrres$method) |
32 | | - { |
33 | | - temp <- mr_egger_regression_bootstrap(d$beta.exposure, d$beta.outcome, d$se.exposure, d$se.outcome, default_parameters()) |
34 | | - mrres$a[mrres$method == "MR Egger (bootstrap)"] <- temp$b_i |
35 | | - } |
| 36 | + if ("MR Egger (bootstrap)" %in% mrres$method) { |
| 37 | + temp <- mr_egger_regression_bootstrap( |
| 38 | + d$beta.exposure, |
| 39 | + d$beta.outcome, |
| 40 | + d$se.exposure, |
| 41 | + d$se.outcome, |
| 42 | + default_parameters() |
| 43 | + ) |
| 44 | + mrres$a[mrres$method == "MR Egger (bootstrap)"] <- temp$b_i |
| 45 | + } |
36 | 46 |
|
37 | | - ggplot2::ggplot(data=d, ggplot2::aes(x=beta.exposure, y=beta.outcome)) + |
38 | | - ggplot2::geom_errorbar(ggplot2::aes(ymin=beta.outcome-se.outcome, ymax=beta.outcome+se.outcome), colour="grey", width=0) + |
39 | | - ggplot2::geom_errorbarh(ggplot2::aes(xmin=beta.exposure-se.exposure, xmax=beta.exposure+se.exposure), colour="grey", height=0) + |
40 | | - ggplot2::geom_point() + |
41 | | - ggplot2::geom_abline(data=mrres, ggplot2::aes(intercept=a, slope=b, colour=method), show.legend=TRUE) + |
42 | | - ggplot2::scale_colour_manual(values=c("#a6cee3", "#1f78b4", "#b2df8a", "#33a02c", "#fb9a99", "#e31a1c", "#fdbf6f", "#ff7f00", "#cab2d6", "#6a3d9a", "#ffff99", "#b15928")) + |
43 | | - ggplot2::labs(colour="MR Test", x=paste("SNP effect on", d$exposure[1]), y=paste("SNP effect on", d$outcome[1])) + |
44 | | - ggplot2::theme(legend.position="top", legend.direction="vertical") + |
45 | | - ggplot2::guides(colour=ggplot2::guide_legend(ncol=2)) |
46 | | - }) |
47 | | - mrres |
| 47 | + if ("MR GRIP" %in% mrres$method) { |
| 48 | + temp <- mr_grip( |
| 49 | + d$beta.exposure, |
| 50 | + d$beta.outcome, |
| 51 | + d$se.exposure, |
| 52 | + d$se.outcome, |
| 53 | + default_parameters() |
| 54 | + ) |
| 55 | + # keep intercept at 0 because plot on gd versus gp axes |
| 56 | + # mrres$a[mrres$method == "MR GRIP"] <- temp$b_i |
| 57 | + msgtxt <- paste0("Strictly, it is only valid to view the MR-GRIP estimate on the standard MR scatter plot axes when the intercept is zero. The estimated intercept for this model is: ", signif(temp$b_i)) |
| 58 | + message(msgtxt) |
| 59 | + } |
| 60 | + |
| 61 | + ggplot2::ggplot( |
| 62 | + data = d, |
| 63 | + ggplot2::aes(x = beta.exposure, y = beta.outcome) |
| 64 | + ) + |
| 65 | + ggplot2::geom_errorbar( |
| 66 | + ggplot2::aes( |
| 67 | + ymin = beta.outcome - se.outcome, |
| 68 | + ymax = beta.outcome + se.outcome |
| 69 | + ), |
| 70 | + colour = "grey", |
| 71 | + width = 0 |
| 72 | + ) + |
| 73 | + ggplot2::geom_errorbarh( |
| 74 | + ggplot2::aes( |
| 75 | + xmin = beta.exposure - se.exposure, |
| 76 | + xmax = beta.exposure + se.exposure |
| 77 | + ), |
| 78 | + colour = "grey", |
| 79 | + height = 0 |
| 80 | + ) + |
| 81 | + ggplot2::geom_point() + |
| 82 | + ggplot2::geom_abline( |
| 83 | + data = mrres, |
| 84 | + ggplot2::aes(intercept = a, slope = b, colour = method), |
| 85 | + show.legend = TRUE |
| 86 | + ) + |
| 87 | + ggplot2::scale_colour_manual( |
| 88 | + values = c( |
| 89 | + "#a6cee3", |
| 90 | + "#1f78b4", |
| 91 | + "#b2df8a", |
| 92 | + "#33a02c", |
| 93 | + "#fb9a99", |
| 94 | + "#e31a1c", |
| 95 | + "#fdbf6f", |
| 96 | + "#ff7f00", |
| 97 | + "#cab2d6", |
| 98 | + "#6a3d9a", |
| 99 | + "#ffff99", |
| 100 | + "#b15928" |
| 101 | + ) |
| 102 | + ) + |
| 103 | + ggplot2::labs( |
| 104 | + colour = "MR Estimate", |
| 105 | + x = paste("SNP effect on", d$exposure[1]), |
| 106 | + y = paste("SNP effect on", d$outcome[1]) |
| 107 | + ) + |
| 108 | + ggplot2::theme(legend.position = "top", legend.direction = "vertical") + |
| 109 | + ggplot2::guides(colour = ggplot2::guide_legend(ncol = 2)) |
| 110 | + }) |
| 111 | + mrres |
48 | 112 | } |
49 | 113 |
|
50 | 114 |
|
51 | | -blank_plot <- function(message) |
52 | | -{ |
53 | | - ggplot2::ggplot(data.frame(a=0,b=0,n=message)) + |
54 | | - ggplot2::geom_text(ggplot2::aes(x=a,y=b,label=n)) + |
55 | | - ggplot2::labs(x=NULL,y=NULL) + |
56 | | - ggplot2::theme(axis.text=ggplot2::element_blank(), axis.ticks=ggplot2::element_blank()) |
| 115 | +blank_plot <- function(message) { |
| 116 | + ggplot2::ggplot(data.frame(a = 0, b = 0, n = message)) + |
| 117 | + ggplot2::geom_text(ggplot2::aes(x = a, y = b, label = n)) + |
| 118 | + ggplot2::labs(x = NULL, y = NULL) + |
| 119 | + ggplot2::theme( |
| 120 | + axis.text = ggplot2::element_blank(), |
| 121 | + axis.ticks = ggplot2::element_blank() |
| 122 | + ) |
57 | 123 | } |
0 commit comments