brielin/WWER
Welch Weighted Egger Regression

Package index
Search the brielin/WWER package
Vignettes
Functions
31
Source code
4
Man pages
15
Home
/
GitHub
/
brielin/WWER
/

In brielin/WWER: Welch Weighted Egger Regression 

library(tidyverse)
library(foreach)
library(doMC)
library(scales)
library(egg)
library(MendelianRandomization)
library(cause)
library(MRMix)
library(MRPRESSO)
library(rslurm)
library(WWER)

# Some helper functions to assist in simualting and generating plots.

# Note that this closes over a globally defined n_cores variable.
run_sim_niter <- function(name, n_iter, Na, Nb, M, ha, hb, hu, nu, eta, gamma, delta, p, q, r, s, af_ref_file = NULL){
  doMC::registerDoMC(cores = n_cores)
  print(c(n_cores, n_iter, Na, Nb, M, ha, hb, hu, nu, eta, gamma, delta, p, q, r, s))
  res_tibble = foreach(iter = 1:n_iter, .combine = bind_rows) %dopar% {
    run_sim(name, iter, Na, Nb, M, ha, hb, hu, nu, eta, gamma, delta, p, q, r, s, af_ref_file)
  }
  return(res_tibble)
}


se <- function(x, na.rm = FALSE){
  sqrt(var(x, na.rm = na.rm)/length(x))
}


make_summary_tibble <- function(plot_tibble){
  summary_tibble <- plot_tibble %>%
    dplyr::group_by(name, method, mr_method, p_thresh, Na, Nb, M, ha, hb, hu, nu, eta, gamma, delta, q, r, s, Rab, Rba) %>% 
    mutate(stat_ab = mean(pab < 0.05, na.rm=T), stat_ba = mean(pba < 0.05, na.rm=T), stat_both = mean(pab < 0.05 & pba < 0.05, na.rm=T),
           se_stat_ab = se(pab < 0.05, na.rm=T), se_stat_ba = se(pba < 0.05, na.rm=T), se_stat_both = se(pab < 0.05 & pba < 0.05, na.rm=T),
           mae_ab = mean(abs(Rab - Rab_hat), na.rm=T), mae_ba = mean(abs(Rba - Rba_hat), na.rm=T),
           se_mae_ab = se(abs(Rab - Rab_hat), na.rm=T), se_mae_ba = se(abs(Rba - Rba_hat), na.rm=T),
           Rab_hat = mean(Rab_hat, na.rm=T), Rba_hat = mean(Rba_hat, na.rm=T),
           Sab_hat = mean(Sab_hat, na.rm=T), Sba_hat = mean(Sba_hat, na.rm=T),
           runtime = mean(runtime, na.rm=T), rg = mean(rg)) %>%
    ungroup() %>%
    distinct(name, method, mr_method, p_thresh, Na, Nb, M, ha, hb, hu, rg, nu, eta, gamma, delta, q, r, s, Rab, Rba,
             stat_ab, stat_ba, stat_both, se_stat_ab, se_stat_ba, se_stat_both, mae_ab, mae_ba, se_mae_ab, se_mae_ba,
             Rab_hat, Rba_hat, Sab_hat, Sba_hat, runtime) %>%
    mutate(method = if_else(method == "egger_wf_5e-8", "aaa_egger_wf_5e-8", method))
  return(summary_tibble)
}


reverselog_trans <- function(base = exp(1)) {
  trans <- function(x) -log(x, base)
  inv <- function(x) base^(-x)
  trans_new(paste0("reverselog-", format(base)), trans, inv,
            log_breaks(base = base), domain = c(1e-100, Inf))
}

Simulations under the one-way alt with no confounding effect.

n_iter = 100
af_ref_file <- "/gpfs/commons/projects/UKBB/sumstats/ldsc_results/eur_w_ld_chr/19.l2.ldscore.gz"

h = 0.25
M = 500000
gamma_alt = c(rep(sqrt(0.01), 3), rep(sqrt(0.025), 3), rep(sqrt(0.05), 3), rep(sqrt(0.1), 3), rep(sqrt(0.2), 3))
Na_alt = rep(c(100000, 25000, 100000), 5)
Nb_alt = rep(c(100000, 100000, 25000), 5)
M_l = 500
M_h = 2000

M_ll = 2*M_l
p_ll = M_ll/M
r_ll = M_l/M_ll
s_ll = M_l/M_ll
alt_ll <- tibble(
  name = "PolyG: low A, low B.",
  n_iter = n_iter, Na = Na_alt, Nb = Nb_alt, M = M,
  ha = h, hb = h, hu = h, nu = 0, eta = 0, gamma = gamma_alt, delta = 0,
  p = p_ll, q = 0, r = r_ll, s = s_ll, af_ref_file = af_ref_file
)

M_lh = M_l + M_h
p_lh = M_lh/M
r_lh = M_l/M_lh
s_lh = M_h/M_lh
alt_lh <- tibble(
  name = "PolyG: low A, high B",
  n_iter = n_iter, Na = Na_alt, Nb = Nb_alt, M = M,
  ha = h, hb = h, hu = h, nu = 0, eta = 0, gamma = gamma_alt, delta = 0,
  p = p_lh, q = 0, r = r_lh, s = s_lh, af_ref_file = af_ref_file
)

M_hl = M_l + M_h
p_hl = M_hl/M
r_hl = M_h/M_hl
s_hl = M_l/M_hl
alt_hl <- tibble(
  name = "PolyG: high A, low B",
  n_iter = n_iter, Na = Na_alt, Nb = Nb_alt, M = M,
  ha = h, hb = h, hu = h, nu = 0, eta = 0, gamma = gamma_alt, delta = 0,
  p = p_hl, q = 0, r = r_hl, s = s_hl, af_ref_file = af_ref_file
)

M_hh = 2*M_h
p_hh = M_hh/M
r_hh = M_h/M_hh
s_hh = M_h/M_hh
alt_hh <- tibble(
  name = "PolyG: high A, high B",
  n_iter = n_iter, Na = Na_alt, Nb = Nb_alt, M = M,
  ha = h, hb = h, hu = h, nu = 0, eta = 0, gamma = gamma_alt, delta = 0,
  p = p_hh, q = 0, r = r_hh, s = s_hh, af_ref_file = af_ref_file
)

sim_settings_alt_np <- dplyr::bind_rows(alt_ll, alt_lh, alt_hl, alt_hh)

n_cores = 8
sjob_alt_np <- rslurm::slurm_apply(
  run_sim_niter, sim_settings_alt_np, jobname = "bimr_alt_np", nodes = nrow(sim_settings_alt_np),
  cpus_per_node = 1, slurm_options = list(time = "36:00:00", mem = "40G"),
  add_objects = c("n_cores"),
  submit = F)

Simulations under the one-way alt with a confounding effect.

h = 0.25
Na = 100000
Nb = 100000
nu_alt_wp = c(sqrt(0.05), sqrt(0.05), sqrt(0.2), sqrt(0.2))
eta_alt_wp = c(sqrt(0.05), sqrt(0.2), sqrt(0.05), sqrt(0.2))
gamma_alt_wp = sqrt(0.025)

M_l = 500
M_h = 2000

M_new = 500000
M_lll = 3*M_l
p_lll = M_lll/M_new
q_lll = M_l/M_lll
r_lll = M_l/M_lll
s_lll = M_l/M_lll
alt_wp_lll <- tibble(
  name = "PolyG: low A, low B, low U",
  n_iter = n_iter, Na = Na, Nb = Nb, M = M,
  ha = h, hb = h, hu = h, nu = nu_alt_wp, eta = eta_alt_wp, gamma = gamma_alt_wp, delta = 0,
  p = p_lll, q = q_lll, r = r_lll, s = s_lll, af_ref_file = af_ref_file
)

M_hll = 1*M_h + 2*M_l
p_hll = M_hll/M_new
q_hll = M_h/M_hll
r_hll = M_l/M_hll
s_hll = M_l/M_hll
alt_wp_hll <- tibble(
  name = "PolyG: low A, low B, high U",
  n_iter = n_iter, Na = Na, Nb = Nb, M = M,
  ha = h, hb = h, hu = h, nu = nu_alt_wp, eta = eta_alt_wp, gamma = gamma_alt_wp, delta = 0,
  p = p_hll, q = q_hll, r = r_hll, s = s_hll, af_ref_file = af_ref_file
)

M_hlh = 2*M_h + 1*M_l
p_hlh = M_hlh/M_new
q_hlh = M_h/M_hll
r_hlh = M_l/M_hll
s_hlh = M_h/M_hll
alt_wp_hlh <- tibble(
  name = "PolyG: high A, low B, high U",
  n_iter = n_iter, Na = Na, Nb = Nb, M = M,
  ha = h, hb = h, hu = h, nu = nu_alt_wp, eta = eta_alt_wp, gamma = gamma_alt_wp, delta = 0,
  p = p_hlh, q = q_hlh, r = r_hlh, s = s_hlh, af_ref_file = af_ref_file
)

M_llh = 1*M_h + 2*M_l
p_llh = M_llh/M_new
q_llh = M_l/M_llh
r_llh = M_l/M_llh
s_llh = M_h/M_llh
alt_wp_llh <- tibble(
  name = "PolyG: low A, high B, low U",
  n_iter = n_iter, Na = Na, Nb = Nb, M = M,
  ha = h, hb = h, hu = h, nu = nu_alt_wp, eta = eta_alt_wp, gamma = gamma_alt_wp, delta = 0,
  p = p_llh, q = q_llh, r = r_llh, s = s_llh, af_ref_file = af_ref_file
)

M_lhh = 2*M_h + 1*M_l
p_lhh = M_lhh/M_new
q_lhh = M_l/M_lhh
r_lhh = M_h/M_lhh
s_lhh = M_h/M_lhh
alt_wp_lhh <- tibble(
  name = "PolyG: high A, high B, low U",
  n_iter = n_iter, Na = Na, Nb = Nb, M = M,
  ha = h, hb = h, hu = h, nu = nu_alt_wp, eta = eta_alt_wp, gamma = gamma_alt_wp, delta = 0,
  p = p_lhh, q = q_lhh, r = r_lhh, s = s_lhh, af_ref_file = af_ref_file
)

M_hhh = 3*M_h
p_hhh = M_hhh/M_new
q_hhh = M_h/M_hhh
r_hhh = M_h/M_hhh
s_hhh = M_h/M_hhh
alt_wp_hhh <- tibble(
  name = "PolyG: high A, high B, high U",
  n_iter = n_iter, Na = Na, Nb = Nb, M = M,
  ha = h, hb = h, hu = h, nu = nu_alt_wp, eta = eta_alt_wp, gamma = gamma_alt_wp, delta = 0,
  p = p_hhh, q = q_hhh, r = r_hhh, s = s_hhh, af_ref_file = af_ref_file
)

sim_settings_alt_wp <- dplyr::bind_rows(alt_wp_lll, alt_wp_hll, alt_wp_hlh, alt_wp_llh, alt_wp_lhh, alt_wp_hhh)

n_cores = 8
sjob_alt_wp <- rslurm::slurm_apply(
  run_sim_niter, sim_settings_alt_wp, jobname = "bimr_alt_wp", nodes = nrow(sim_settings_alt_wp),
  cpus_per_node = 1, slurm_options = list(time = "36:00:00", mem = "40G"),
  add_objects = c("n_cores"),
  submit = F)

Simulations under the bi-directional alt.

n_iter = 100
af_ref_file <- "/gpfs/commons/projects/UKBB/sumstats/ldsc_results/eur_w_ld_chr/19.l2.ldscore.gz"

h = 0.25
M = 500000
# Originally I had negative delta, but results are symmetric in delta.
gamma_alt = c(rep(-sqrt(0.1), 6), rep(-sqrt(0.03), 6), rep(-sqrt(0.01), 6), rep(sqrt(0.01), 6), rep(sqrt(0.03), 6), rep(sqrt(0.1), 6))
delta_alt = rep(c(-sqrt(0.1), -sqrt(0.03), -sqrt(0.01), sqrt(0.01), sqrt(0.03), sqrt(0.1)), 6)
# delta_alt = rep(c(sqrt(0.01), sqrt(0.03), sqrt(0.1)), 3)
# gamma_alt = c(rep(-sqrt(0.1), 3), rep(-sqrt(0.03), 3), rep(-sqrt(0.01), 3), rep(sqrt(0.01), 3), rep(sqrt(0.03), 3), rep(sqrt(0.1), 3))
M_l = 500
M_h = 2000

M_hh = 2*M_h
p_hh = M_hh/M
r_hh = M_h/M_hh
s_hh = M_h/M_hh
bialt_hh_eq <- tibble(
  name = "PolyG: high A, high B; Na = Nb",
  n_iter = n_iter, Na = 100000, Nb = 100000, M = M,
  ha = h, hb = h, hu = h, nu = 0, eta = 0, gamma = gamma_alt, delta = delta_alt,
  p = p_hh, q = 0, r = r_hh, s = s_hh, af_ref_file = af_ref_file
)

M_hh = 2*M_h
p_hh = M_hh/M
r_hh = M_h/M_hh
s_hh = M_h/M_hh
bialt_hh_agb <- tibble(
  name = "PolyG: high A, high B; Na > Nb",
  n_iter = n_iter, Na = 100000, Nb = 25000, M = M,
  ha = h, hb = h, hu = h, nu = 0, eta = 0, gamma = gamma_alt, delta = delta_alt,
  p = p_hh, q = 0, r = r_hh, s = s_hh, af_ref_file = af_ref_file
)

M_hl = M_l + M_h
p_hl = M_hl/M
r_hl = M_h/M_hl
s_hl = M_l/M_hl
bialt_hl_eq <- tibble(
  name = "PolyG: high A, low B; Na = Nb",
  n_iter = n_iter, Na = 100000, Nb = 100000, M = M,
  ha = h, hb = h, hu = h, nu = 0, eta = 0, gamma = gamma_alt, delta = delta_alt,
  p = p_hl, q = 0, r = r_hl, s = s_hl, af_ref_file = af_ref_file
)

M_hl = M_l + M_h
p_hl = M_hl/M
r_hl = M_h/M_hl
s_hl = M_l/M_hl
bialt_hl_agb <- tibble(
  name = "PolyG: high A, low B; Na > Nb",
  n_iter = n_iter, Na = 100000, Nb = 25000, M = M,
  ha = h, hb = h, hu = h, nu = 0, eta = 0, gamma = gamma_alt, delta = delta_alt,
  p = p_hl, q = 0, r = r_hl, s = s_hl, af_ref_file = af_ref_file
)

M_hl = M_l + M_h
p_hl = M_hl/M
r_hl = M_h/M_hl
s_hl = M_l/M_hl
bialt_hl_bga <- tibble(
  name = "PolyG: high A, low B; Na < Nb",
  n_iter = n_iter, Na = 25000, Nb = 100000, M = M,
  ha = h, hb = h, hu = h, nu = 0, eta = 0, gamma = gamma_alt, delta = delta_alt,
  p = p_hl, q = 0, r = r_hl, s = s_hl, af_ref_file = af_ref_file
)

sim_settings_bi_alt <- dplyr::bind_rows(bialt_hh_eq, bialt_hh_agb, bialt_hl_eq, bialt_hl_agb, bialt_hl_bga)

n_cores = 8
sjob_bi_alt <- rslurm::slurm_apply(
  run_sim_niter, sim_settings_bi_alt, jobname = "bimr_bi_alt", nodes = nrow(sim_settings_bi_alt),
  cpus_per_node = 1, slurm_options = list(time = "36:00:00", mem = "40G"),
  add_objects = c("n_cores"),
  submit = F)

sjob_alt_np <- rslurm::slurm_job(jobname = "bimr_alt_np", nodes = nrow(sim_settings_alt_np))
sjob_alt_wp <- rslurm::slurm_job(jobname = "bimr_alt_wp", nodes = nrow(sim_settings_alt_wp))
sjob_bi_alt <- rslurm::slurm_job(jobname = "bimr_bi_alt", nodes = nrow(sim_settings_bi_alt))

alt_np_results <- rslurm::get_slurm_out(sjob_alt_np, outtype = "table")
alt_wp_results <- rslurm::get_slurm_out(sjob_alt_wp, outtype = "table")
bi_alt_results <- rslurm::get_slurm_out(sjob_bi_alt, outtype = "table")

alt_np_results_summary <- make_summary_tibble(alt_np_results)
alt_wp_results_summary <- make_summary_tibble(alt_wp_results)
bi_alt_results_summary <- make_summary_tibble(bi_alt_results)

# Figure 3 for np alt and bi-alt
theme_set(theme_classic())

# show_methods = c("egger_5e-8", "egger_wf_5e-8", "egger_s_5e-8", "cause_1e-3", "ivw_5e-8", "mr_mix_5e-8", "mr_presso_5e-8")
show_methods = c("egger_5e-8", "aaa_egger_wf_5e-8", "egger_s_5e-8", "cause_1e-3", "ivw_5e-8", "mr_mix_5e-8", "mr_presso_5e-8", "mbe_5e-8")
show_labels = c("WWER", "CAUSE", "Egger", "Steiger", "IVW", "MBE", "MR Mix", "MR PRESSO")

g1 <- alt_np_results_summary %>%
  filter(method %in% show_methods) %>%
  filter(Na == 100000, Nb == 100000) %>%
  ggplot() +
  geom_line(aes(x=gamma, y=stat_ab, color=method)) +
  geom_errorbar(aes(ymin = stat_ab - 1.96*se_stat_ab, ymax = stat_ab + 1.96*se_stat_ab, x = gamma), width = 0.025) +
  facet_wrap(vars(name), nrow=1) +
  scale_color_brewer(labels = show_labels, palette = "Dark2") +
  labs(x = NULL, title = "A) One way alt, equal sample size", y = NULL) + guides(color = FALSE)

g2 <- alt_np_results_summary %>%
  filter(method %in% show_methods) %>%
  filter(Na == 25000, Nb == 100000) %>%
  ggplot() +
  geom_line(aes(x=gamma, y=stat_ab, color=method)) +
  geom_errorbar(aes(ymin = stat_ab - 1.96*se_stat_ab, ymax = stat_ab + 1.96*se_stat_ab, x = gamma), width = 0.025) +
  facet_wrap(vars(name), nrow=1) +
  scale_color_brewer(labels = show_labels, palette = "Dark2") +
  labs(x = NULL, title = "B) One way alt, Na < Nb", y = NULL) + guides(color = FALSE)

g3 <- alt_np_results_summary %>%
  filter(method %in% show_methods) %>%
  filter(Na == 100000, Nb == 25000) %>%
  ggplot() +
  geom_line(aes(x=gamma, y=stat_ab, color=method)) +
  geom_errorbar(aes(ymin = stat_ab - 1.96*se_stat_ab, ymax = stat_ab + 1.96*se_stat_ab, x = gamma), width = 0.025) +
  facet_wrap(vars(name), nrow=1) +
  scale_color_brewer(labels = show_labels, palette = "Dark2") +
  labs(x = NULL, title = "C) One way alt, Na > Nb", y = NULL, color = "Method")

g4 <- bi_alt_results_summary %>%
  filter(delta == sqrt(0.01)) %>%
  filter(method %in% show_methods) %>%
  ggplot() +
  geom_line(aes(x=gamma, y=stat_both, color=method)) +
  geom_errorbar(aes(ymin = stat_both - 1.96*se_stat_both, ymax = stat_both + 1.96*se_stat_both, x = gamma), width = 0.025) +
  facet_wrap(vars(name), nrow=1) +
  scale_color_brewer(labels = show_labels, palette = "Dark2") +
  labs(x = NULL, title = "D) Bi-directional alt, delta = 0.1", y = NULL) + guides(color = FALSE)


g5 <- bi_alt_results_summary %>%
  filter(delta == sqrt(0.03)) %>%
  filter(method %in% show_methods) %>%
  ggplot() +
  geom_line(aes(x=gamma, y=stat_both, color=method)) +
  geom_errorbar(aes(ymin = stat_both - 1.96*se_stat_both, ymax = stat_both + 1.96*se_stat_both, x = gamma), width = 0.025) +
  facet_wrap(vars(name), nrow=1) +
  scale_color_brewer(labels = show_labels, palette = "Dark2") +
  labs(x = NULL, title = "E) Bi-directional alt, delta = 0.17", y = NULL) + guides(color = FALSE)


g6 <- bi_alt_results_summary %>%
  filter(delta == sqrt(0.1)) %>%
  filter(method %in% show_methods) %>%
  ggplot() +
  geom_line(aes(x=gamma, y=stat_both, color=method)) +
  geom_errorbar(aes(ymin = stat_both - 1.96*se_stat_both, ymax = stat_both + 1.96*se_stat_both, x = gamma), width = 0.025) +
  facet_wrap(vars(name), nrow=1) +
  scale_color_brewer(labels = show_labels, palette = "Dark2") +
  labs(x = NULL, title = "F) Bi-directional alt, delta = 0.32", y = NULL) + guides(color = FALSE)

plot = ggarrange(g1, g2, g3, g4, g5, g6, nrow=6, left = "Power", bottom = "Gamma")
ggsave("../Figures/wwer_fig3.pdf", device = "pdf", plot = plot, dpi = 300)

# Supp figure with reverse direction FPR
g1 <- alt_np_results_summary %>%
  filter(method %in% show_methods) %>%
  filter(Na == 100000, Nb == 25000) %>%
  ggplot() +
  geom_line(aes(x=gamma, y=stat_ba, color=method)) +
  geom_errorbar(aes(ymin = stat_ba - 1.96*se_stat_ba, ymax = stat_ba + 1.96*se_stat_ba, x = gamma), width = 0.025) +
  facet_wrap(vars(name), nrow=1) +
  ylim(0, 1) +
  scale_color_brewer(labels =show_labels, palette = "Dark2") +
  labs(x = NULL, title = "C) One way alt reverse positivity, study A higher power", y = NULL) + guides(color = FALSE)

g2 <- alt_np_results_summary %>%
  filter(method %in% show_methods) %>%
  filter(Na == 25000, Nb == 100000) %>%
  ggplot() +
  geom_line(aes(x=gamma, y=stat_ba, color=method)) +
  geom_errorbar(aes(ymin = stat_ba - 1.96*se_stat_ba, ymax = stat_ba + 1.96*se_stat_ba, x = gamma), width = 0.025) +
  facet_wrap(vars(name), nrow=1) +
  scale_color_brewer(labels =show_labels, palette = "Dark2") +
  labs(x = NULL, title = "B) One way alt reverse positivity, study B higher power", y = NULL, color = "Method")

g3 <- alt_np_results_summary %>%
  filter(method %in% show_methods) %>%
  filter(Na == 100000, Nb == 100000) %>%
  ggplot() +
  geom_line(aes(x=gamma, y=stat_ba, color=method)) +
  geom_errorbar(aes(ymin = stat_ba - 1.96*se_stat_ba, ymax = stat_ba + 1.96*se_stat_ba, x = gamma), width = 0.025) +
  facet_wrap(vars(name), nrow=1) +
  scale_color_brewer(labels =show_labels, palette = "Dark2") +
  labs(x = NULL, title = "A) One way alt reverse positivity, equal sample sizes", y = NULL) + guides(color = FALSE)

ggarrange(g3, g2, g1, nrow=3, left = "FPR, B to A direction", bottom = "Gamma, A to B direction")

g1 <- alt_wp_results_summary %>%
  filter(nu == sqrt(0.05), eta == sqrt(0.05)) %>%
  filter(method %in% show_methods) %>%
  ggplot() +
  geom_col(aes(x=method, y=stat_ab, fill=method), position = "dodge") + 
  geom_errorbar(aes(ymin = stat_ab - 1.96*se_stat_ab, ymax = stat_ab + 1.96*se_stat_ab, x=method), width = 0.15) +
  scale_fill_brewer(labels =show_labels, palette = "Dark2") +
  ylim(0, 1) +
  theme(axis.text.x=element_blank(), axis.ticks.x=element_blank()) +
  facet_wrap(vars(name), nrow=1) +
  labs(x = NULL, title = "A) One way alt with equal weak pleiotropy", y = NULL) + guides(fill = FALSE)

g2 <- alt_wp_results_summary %>%
  filter(nu == sqrt(0.2), eta == sqrt(0.2)) %>%
  filter(method %in% show_methods) %>%
  ggplot() +
  geom_col(aes(x=method, y=stat_ab, fill=method), position = "dodge") + 
  geom_errorbar(aes(ymin = stat_ab - 1.96*se_stat_ab, ymax = stat_ab + 1.96*se_stat_ab, x=method), width = 0.15) +
  scale_fill_brewer(labels =show_labels, palette = "Dark2") +
  ylim(0, 1) +
  theme(axis.text.x=element_blank(), axis.ticks.x=element_blank()) +
  facet_wrap(vars(name), nrow=1) +
  labs(x = NULL, title = "B) One way alt with equal strong pleiotropy", y = NULL, fill = "Method")

g3 <- alt_wp_results_summary %>%
  filter(nu == sqrt(0.05), eta == sqrt(0.2)) %>%
  filter(method %in% show_methods) %>%
  ggplot() +
  geom_col(aes(x=method, y=stat_ab, fill=method), position = "dodge") + 
  geom_errorbar(aes(ymin = stat_ab - 1.96*se_stat_ab, ymax = stat_ab + 1.96*se_stat_ab, x=method), width = 0.15) +
  scale_fill_brewer(labels =show_labels, palette = "Dark2") +
  ylim(0, 1) +
  theme(axis.text.x=element_blank(), axis.ticks.x=element_blank()) +
  facet_wrap(vars(name), nrow=1) +
  labs(x = NULL, title = "C) One way alt with stronger pleiotropic effect on A", y = NULL) + guides(fill = FALSE)

g4 <- alt_wp_results_summary %>%
  filter(nu == sqrt(0.2), eta == sqrt(0.05)) %>%
  filter(method %in% show_methods) %>%
  ggplot() +
  geom_col(aes(x=method, y=stat_ab, fill=method), position = "dodge") + 
  geom_errorbar(aes(ymin = stat_ab - 1.96*se_stat_ab, ymax = stat_ab + 1.96*se_stat_ab, x=method), width = 0.15) +
  scale_fill_brewer(labels =show_labels, palette = "Dark2") +
  ylim(0, 1) +
  theme(axis.text.x=element_blank(), axis.ticks.x=element_blank()) +
  facet_wrap(vars(name), nrow=1) +
  labs(x = NULL, title = "D) One way alt with stronger pleiotropic effect on B", y = NULL) + guides(fill = FALSE)

ggarrange(g1, g2, g3, g4, nrow=4, left = "Power, A to B direction", bottom = "Method")

g1 <- alt_wp_results_summary %>%
  filter(nu == sqrt(0.05), eta == sqrt(0.05)) %>%
  filter(method %in% show_methods) %>%
  ggplot() +
  geom_col(aes(x=method, y=stat_ba, fill=method), position = "dodge") + 
  geom_errorbar(aes(ymin = stat_ba - 1.96*se_stat_ba, ymax = stat_ba + 1.96*se_stat_ba, x=method), width = 0.15) +
  scale_fill_brewer(labels =show_labels, palette = "Dark2") +
  ylim(0, 1) +
  theme(axis.text.x=element_blank(), axis.ticks.x=element_blank()) +
  facet_wrap(vars(name), nrow=1) +
  labs(x = NULL, title = "A) One way alt with equal weak pleiotropy", y = NULL) + guides(fill = FALSE)

g2 <- alt_wp_results_summary %>%
  filter(nu == sqrt(0.2), eta == sqrt(0.2)) %>%
  filter(method %in% show_methods) %>%
  ggplot() +
  geom_col(aes(x=method, y=stat_ba, fill=method), position = "dodge") + 
  geom_errorbar(aes(ymin = stat_ba - 1.96*se_stat_ba, ymax = stat_ba + 1.96*se_stat_ba, x=method), width = 0.15) +
  scale_fill_brewer(labels =show_labels, palette = "Dark2") +
  ylim(0, 1) +
  theme(axis.text.x=element_blank(), axis.ticks.x=element_blank()) +
  facet_wrap(vars(name), nrow=1) +
  labs(x = NULL, title = "B) One way alt with equal strong pleiotropy", y = NULL, fill = "Method")

g3 <- alt_wp_results_summary %>%
  filter(nu == sqrt(0.05), eta == sqrt(0.2)) %>%
  filter(method %in% show_methods) %>%
  ggplot() +
  geom_col(aes(x=method, y=stat_ba, fill=method), position = "dodge") + 
  geom_errorbar(aes(ymin = stat_ba - 1.96*se_stat_ba, ymax = stat_ba + 1.96*se_stat_ba, x=method), width = 0.15) +
  scale_fill_brewer(labels =show_labels, palette = "Dark2") +
  ylim(0, 1) +
  theme(axis.text.x=element_blank(), axis.ticks.x=element_blank()) +
  facet_wrap(vars(name), nrow=1) +
  labs(x = NULL, title = "C) One way alt with stronger pleiotropic effect on A", y = NULL) + guides(fill = FALSE)

g4 <- alt_wp_results_summary %>%
  filter(nu == sqrt(0.2), eta == sqrt(0.05)) %>%
  filter(method %in% show_methods) %>%
  ggplot() +
  geom_col(aes(x=method, y=stat_ba, fill=method), position = "dodge") + 
  geom_errorbar(aes(ymin = stat_ba - 1.96*se_stat_ba, ymax = stat_ba + 1.96*se_stat_ba, x=method), width = 0.15) +
  scale_fill_brewer(labels =show_labels, palette = "Dark2") +
  ylim(0, 1) +
  theme(axis.text.x=element_blank(), axis.ticks.x=element_blank()) +
  facet_wrap(vars(name), nrow=1) +
  labs(x = NULL, title = "D) One way alt with stronger pleiotropic effect on B", y = NULL) + guides(fill = FALSE)

ggarrange(g1, g2, g3, g4, nrow=4, left = "FPR, B to A direction", bottom = "Method")

# Table time
res_pow_np_ab <- alt_np_results_summary %>% dplyr::filter(p_thresh != "5e-6", mr_method != "egger_wa") %>%
  dplyr::mutate(POW = paste0(stat_ab, " (", round(se_stat_ab, digits=2), ")"), study = paste(name, round(Na, 2), round(Nb, 2), round(gamma, 2))) %>%
  dplyr::select(study, mr_method, POW) %>% tidyr::pivot_wider(names_from = mr_method, values_from = POW)

res_fpr_np_ba <- alt_np_results_summary %>% dplyr::filter(p_thresh != "5e-6", mr_method != "egger_wa") %>%
  dplyr::mutate(FPR = paste0(stat_ba, " (", round(se_stat_ba, digits=2), ")"), study = paste(name, round(Na, 2), round(Nb, 2), round(gamma, 2))) %>%
  dplyr::select(study, mr_method, FPR) %>% tidyr::pivot_wider(names_from = mr_method, values_from = FPR)

res_pow_wp_ab <- alt_wp_results_summary %>% dplyr::filter(p_thresh != "5e-6", mr_method != "egger_wa") %>%
  dplyr::mutate(POW = paste0(stat_ab, " (", round(se_stat_ab, digits=2), ")"), study = paste(name, round(nu, 2), round(eta, 2))) %>%
  dplyr::select(study, mr_method, POW) %>% tidyr::pivot_wider(names_from = mr_method, values_from = POW)

res_fpr_wp_ba <- alt_wp_results_summary %>% dplyr::filter(p_thresh != "5e-6", mr_method != "egger_wa") %>%
  dplyr::mutate(FPR = paste0(stat_ba, " (", round(se_stat_ba, digits=2), ")"), study = paste(name, round(nu, 2), round(eta, 2))) %>%
  dplyr::select(study, mr_method, FPR) %>% tidyr::pivot_wider(names_from = mr_method, values_from = FPR)

res_pow_bi_alt <- bi_alt_results_summary %>% dplyr::filter(p_thresh != "5e-6", mr_method != "egger_wa") %>%
  dplyr::mutate(POW = paste0(stat_both, " (", round(se_stat_both, digits=2), ")"), study = paste(name, round(gamma, 2), round(delta, 2))) %>%
  dplyr::select(study, mr_method, POW) %>% tidyr::pivot_wider(names_from = mr_method, values_from = POW)

outfile <- "supp_tables_alt.xlsx"

supp_tables <- openxlsx::createWorkbook()
openxlsx::addWorksheet(supp_tables, "ST10")
openxlsx::addWorksheet(supp_tables, "ST11")
openxlsx::addWorksheet(supp_tables, "ST12")
openxlsx::addWorksheet(supp_tables, "ST13")
openxlsx::addWorksheet(supp_tables, "ST14")
openxlsx::addWorksheet(supp_tables, "ST15")
openxlsx::addWorksheet(supp_tables, "ST16")
openxlsx::addWorksheet(supp_tables, "ST17")
openxlsx::writeData(supp_tables, "ST10", sim_settings_alt_np)
openxlsx::writeData(supp_tables, "ST11", sim_settings_alt_wp)
openxlsx::writeData(supp_tables, "ST12", sim_settings_bi_alt)
openxlsx::writeData(supp_tables, "ST13", res_pow_np_ab)
openxlsx::writeData(supp_tables, "ST14", res_fpr_np_ba)
openxlsx::writeData(supp_tables, "ST15", res_pow_wp_ab)
openxlsx::writeData(supp_tables, "ST16", res_fpr_wp_ba)
openxlsx::writeData(supp_tables, "ST17", res_pow_bi_alt)


openxlsx::saveWorkbook(wb = supp_tables, file = outfile, overwrite = TRUE)


brielin/WWER documentation built on May 22, 2022, 7:28 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close