run_sim <- function(nsnp, n, h2_x, h2_y, h2_u, b_us, b_xs, b_ints, b_uy, sim=1) {
args <- tibble(nsnp, n, h2_x, h2_u, b_us, b_xs, b_ints, b_uy, sim)
gu <- make_geno(n, nsnp, 0.5)
effsu <- choose_effects(nsnp, h2_u)
u <- make_phen(effsu, gu)
uxy <- rnorm(n)
gx <- make_geno(n, nsnp, 0.5)
effsx <- choose_effects(nsnp, h2_x)
x <- make_phen(c(0.2, effsx), cbind(uxy, gx))
test <- rbinom(n, 1, plogis(u * b_us + x * b_xs + u*x * b_ints))
y <- make_phen(c(b_uy, 0.2), cbind(u, uxy))
res <- list()
res[[1]] <- summary(lm(y ~ x))$coef %>% as_tibble() %>% janitor::clean_names() %>% mutate(method="obs", exp="x", out="y", sel="none") %>% slice(2)
res[[2]] <- summary(lm(y[test==1] ~ x[test==1]))$coef %>% as_tibble() %>% janitor::clean_names() %>% mutate(method="obs", exp="x", out="y", sel="x and y") %>% slice(2)
gwas_x <- gwas(x, cbind(gu, gx))
gwas_y <- gwas(y, cbind(gu, gx))
sel <- gwas_x$pval < 5e-8
res[[3]] <- summary(lm(gwas_y$bhat[sel] ~ 0 + gwas_x$bhat[sel], weight=1/gwas_y$se[sel]^2 ))$coef %>% as_tibble() %>% janitor::clean_names() %>% mutate(method="mr", exp="x", out="y", sel="none") %>% slice(1)
gwas_x_test <- gwas(x[test==1], cbind(gu, gx)[test==1,])
gwas_y_test <- gwas(y[test==1], cbind(gu, gx)[test==1,])
sel <- gwas_x_test$pval < 5e-8
res[[4]] <- summary(lm(gwas_y_test$bhat[sel] ~ 0 + gwas_x_test$bhat[sel], weight=1/gwas_y_test$se[sel]^2 ))$coef %>% as_tibble() %>% janitor::clean_names() %>% mutate(method="mr", exp="x", out="y", sel="x and y") %>% slice(1)
res[[5]] <- summary(lm(gwas_y$bhat[sel] ~ 0 + gwas_x_test$bhat[sel], weight=1/gwas_y$se[sel]^2 ))$coef %>% as_tibble() %>% janitor::clean_names() %>% mutate(method="mr", exp="x", out="y", sel="x") %>% slice(1)
res[[6]] <- summary(lm(gwas_y_test$bhat[sel] ~ 0 + gwas_x$bhat[sel], weight=1/gwas_y_test$se[sel]^2 ))$coef %>% as_tibble() %>% janitor::clean_names() %>% mutate(method="mr", exp="x", out="y", sel="y") %>% slice(1)
gwas_test <- gwas(test, cbind(gu, gx))
sel <- gwas_x_test$pval < 5e-8 | gwas_test$pval < 5e-8
res[[7]] <- summary(lm(gwas_y_test$bhat[sel] ~ 0 + gwas_x_test$bhat[sel] + gwas_test$bhat[sel], weight=1/gwas_y_test$se[sel]^2))$coef %>% as_tibble() %>% janitor::clean_names() %>% mutate(method="mvmr", exp="x", out="y", sel="x and y") %>% slice(1)
res[[8]] <- summary(lm(gwas_y$bhat[sel] ~ 0 + gwas_x_test$bhat[sel] + gwas_test$bhat[sel], weight=1/gwas_y$se[sel]^2))$coef %>% as_tibble() %>% janitor::clean_names() %>% mutate(method="mvmr", exp="x", out="y", sel="x") %>% slice(1)
res[[9]] <- summary(lm(gwas_y_test$bhat[sel] ~ 0 + gwas_x$bhat[sel] + gwas_test$bhat[sel], weight=1/gwas_y_test$se[sel]^2))$coef %>% as_tibble() %>% janitor::clean_names() %>% mutate(method="mvmr", exp="x", out="y", sel="y") %>% slice(1)
res[[10]] <- summary(lm(gwas_y$bhat[sel] ~ 0 + gwas_x$bhat[sel] + gwas_test$bhat[sel], weight=1/gwas_y$se[sel]^2))$coef %>% as_tibble() %>% janitor::clean_names() %>% mutate(method="mvmr", exp="x", out="y", sel="none") %>% slice(1)
return(bind_rows(res) %>% select(method, exp, out, sel, everything()) %>% bind_cols(args, .))
}
