From 47b00c9283e09820818d1ddac504870e0b3558b9 Mon Sep 17 00:00:00 2001 From: jiayi-s Date: Thu, 30 Nov 2023 09:12:15 -0500 Subject: [PATCH] added use_robust_var_est option in mJAM_build_CS --- R/mJAM_Forward.R | 3 ++- R/mJAM_build_CS.R | 55 ++++++++++++++++++++++++++++------------ man/mJAM_build_CS.Rd | 5 +++- man/mJAM_get_PrM.Rd | 5 +++- man/mJAM_get_PrM_Wald.Rd | 5 +++- 5 files changed, 53 insertions(+), 20 deletions(-) diff --git a/R/mJAM_Forward.R b/R/mJAM_Forward.R index 4d35a4d..f5244eb 100644 --- a/R/mJAM_Forward.R +++ b/R/mJAM_Forward.R @@ -328,7 +328,8 @@ mJAM_Forward <- function(N_GWAS, X_ref, GItGI_curr = GItGI_curr, GIty_curr = GIty_curr, yty_curr = yty_curr, yty_med = yty_med, N_GWAS = N_GWAS, rare_SNPs = rare_SNPs, - Pr_Med_cut = Pr_Med_cut) + Pr_Med_cut = Pr_Med_cut, + use_robust_var_est = use_robust_var_est) ## step 3: prune out CS SNPs and SNPs in LD with index SNP; subset input statistics curr_CS_snp <- filter(curr_CS, CS_in == T) %>% pull(CS_SNP) diff --git a/R/mJAM_build_CS.R b/R/mJAM_build_CS.R index 80a1770..14af0b0 100644 --- a/R/mJAM_build_CS.R +++ b/R/mJAM_build_CS.R @@ -12,6 +12,7 @@ #' @param prev_X_list A numeric vector of the ID(s) of previously selected index SNP(s). #' @param g The pre-specified `g` in `g`-prior formulation. #' @param rare_SNPs A numeric vector of ID(s) for rare SNP(s) which we do not apply weighting. Instead, we use the individual estimate of yty for these SNPs for robustness. +#' @param use_robust_var_est whether to use linear combination of median yty and individual yty. #' #' #' @author Jiayi Shen @@ -28,7 +29,8 @@ #' } #' -mJAM_get_PrM <- function(GItGI,GIty,yty,yty_med,N_GWAS, C_id,prev_X_list = NULL, g = NULL,rare_SNPs = NULL){ +mJAM_get_PrM <- function(GItGI,GIty,yty,yty_med,N_GWAS, C_id,prev_X_list = NULL, g = NULL,rare_SNPs = NULL, + use_robust_var_est = FALSE){ numEthnic <- length(GItGI) if(is.null(g)){g <- sum(N_GWAS)} @@ -58,9 +60,14 @@ mJAM_get_PrM <- function(GItGI,GIty,yty,yty_med,N_GWAS, C_id,prev_X_list = NULL, }else{ w[[i]] <- yty_med[[i]]/(yty_med[[i]]+delta[[i]]) } - yty_C[[i]] <- w[[i]]*yty[[i]][C_id] + (1-w[[i]])*yty_med[[i]] - # yty_C[[i]] <- yty_med[[i]] # (v4.1a) - # yty_C[[i]] <-yty[[i]][C_id] # (v4.1c) + if(use_robust_var_est){ + yty_C[[i]] <- w[[i]]*yty[[i]][C_id] + (1-w[[i]])*yty_med[[i]] + }else{ + yty_C[[i]] <-yty[[i]][C_id] # (v4.1c) + # yty_C[[i]] <- yty_med[[i]] # (v4.1a) + } + + } susie_in_XtX <- Reduce("+", GItGI_C) susie_in_Xty <- Reduce("+", GIty_C) @@ -120,6 +127,8 @@ mJAM_get_PrM <- function(GItGI,GIty,yty,yty_med,N_GWAS, C_id,prev_X_list = NULL, return(list(post_prob=post_prob,R2_est = R2_est, n_miss = length(missing_ethnic_idx))) } + + #' Get Pr(Model) based on Wald-type model probability #' #' @description Also apply weighting to get robust estimates of yty @@ -133,6 +142,7 @@ mJAM_get_PrM <- function(GItGI,GIty,yty,yty_med,N_GWAS, C_id,prev_X_list = NULL, #' @param prev_X_list A numeric vector of the ID(s) of previously selected index SNP(s). #' @param g The pre-specified `g` in `g`-prior formulation. #' @param rare_SNPs A numeric vector of ID(s) for rare SNP(s) which we do not apply weighting. Instead, we use the individual estimate of yty for these SNPs for robustness. +#' @param use_robust_var_est whether to use linear combination of median yty and individual yty. #' #' #' @author Jiayi Shen @@ -143,7 +153,8 @@ mJAM_get_PrM <- function(GItGI,GIty,yty,yty_med,N_GWAS, C_id,prev_X_list = NULL, #' #' -mJAM_get_PrM_Wald <- function(GItGI, GIty, yty, yty_med, N_GWAS, C_id,prev_X_list = NULL, g = NULL, rare_SNPs = NULL){ +mJAM_get_PrM_Wald <- function(GItGI, GIty, yty, yty_med, N_GWAS, C_id,prev_X_list = NULL, g = NULL, rare_SNPs = NULL, + use_robust_var_est = FALSE){ ## C -> Y | previous index SNPs @@ -163,7 +174,12 @@ mJAM_get_PrM_Wald <- function(GItGI, GIty, yty, yty_med, N_GWAS, C_id,prev_X_lis }else{ w[[i]] <- yty_med[[i]]/(yty_med[[i]]+delta[[i]]) } - yty_C[[i]] <- w[[i]]*yty[[i]][C_id] + (1-w[[i]])*yty_med[[i]] + if(use_robust_var_est){ + yty_C[[i]] <- w[[i]]*yty[[i]][C_id] + (1-w[[i]])*yty_med[[i]] + }else{ + yty_C[[i]] <-yty[[i]][C_id] # (v4.1c) + # yty_C[[i]] <- yty_med[[i]] # (v4.1a) + } } multi_GItGI_C <- Reduce("+", GItGI_C) @@ -319,6 +335,7 @@ mJAM_get_PrMed <- function(GItGI, GIty,yty, yty_med, N_GWAS, g = NULL,C_id, X_id #' @param N_GWAS A vector of sample sizes in all original GWAS studies. #' @param rare_SNPs A numeric vector of ID(s) for rare SNP(s) which we do not apply weighting. Instead, we use the individual estimate of yty for these SNPs for robustness. #' @param Pr_Med_cut The cutoff for Pr(Mediation); SNPs with Pr(Mediation) smaller than this cutoff will be assigned a Pr(CS) = 0 and thus not included in the credible set for the current index +#' @param use_robust_var_est whether to use linear combination of median yty and individual yty. #' #' #' @author Jiayi Shen @@ -341,7 +358,8 @@ mJAM_get_PrMed <- function(GItGI, GIty,yty, yty_med, N_GWAS, g = NULL,C_id, X_id #' } #' mJAM_build_CS <- function(X_id, prev_X_list = NULL,All_id, PrCS_weights = "Pr(M_C)", coverage = 0.95, - GItGI_curr, GIty_curr, yty_curr, yty_med, N_GWAS,rare_SNPs = NULL, Pr_Med_cut = 0.1){ + GItGI_curr, GIty_curr, yty_curr, yty_med, N_GWAS,rare_SNPs = NULL, Pr_Med_cut = 0.1, + use_robust_var_est = FALSE){ ############################################################################################### ## --- Initiate all variables @@ -360,7 +378,8 @@ mJAM_build_CS <- function(X_id, prev_X_list = NULL,All_id, PrCS_weights = "Pr(M_ yty = yty_curr, yty_med = yty_med, N_GWAS = N_GWAS, C_id = X_id, prev_X_list = prev_X_list, g = sum(N_GWAS), - rare_SNPs = rare_SNPs) + rare_SNPs = rare_SNPs, + use_robust_var_est = use_robust_var_est) Post_Model_Prob[X_id] <- Post_CS_Prob[X_id] <- sum_Post_CS_Porb <- Index_Model_Prob <- temp_PrM_res$post_prob Post_Model_Prob_Ratio[X_id] <- 1 } @@ -368,7 +387,8 @@ mJAM_build_CS <- function(X_id, prev_X_list = NULL,All_id, PrCS_weights = "Pr(M_ Post_Model_Prob[X_id] <- Post_CS_Prob[X_id] <- sum_Post_CS_Porb <- Index_Model_Prob <- mJAM_get_PrM_Wald(GItGI = GItGI_curr, GIty = GIty_curr, yty = yty_curr, yty_med = yty_med, N_GWAS = N_GWAS, C_id = X_id, rare_SNPs = rare_SNPs, - prev_X_list = prev_X_list, g = sum(N_GWAS)) + prev_X_list = prev_X_list, g = sum(N_GWAS), + use_robust_var_est = use_robust_var_est) Post_Model_Prob_Ratio[X_id] <- 1 } @@ -388,13 +408,15 @@ mJAM_build_CS <- function(X_id, prev_X_list = NULL,All_id, PrCS_weights = "Pr(M_ yty = yty_curr, yty_med = yty_med, N_GWAS = N_GWAS, C_id = C_id, prev_X_list = prev_X_list, g = sum(N_GWAS), - rare_SNPs = rare_SNPs) + rare_SNPs = rare_SNPs, + use_robust_var_est = use_robust_var_est) temp_r2 <- temp_PrM_res$post_prob } if(PrCS_weights == "Pr(Wald)"){ temp_r2 <- mJAM_get_PrM_Wald(GItGI = GItGI_curr, GIty = GIty_curr, N_GWAS = N_GWAS, yty = yty_curr, yty_med = yty_med, rare_SNPs = rare_SNPs, - C_id = C_id, prev_X_list = prev_X_list, g = sum(N_GWAS)) + C_id = C_id, prev_X_list = prev_X_list, g = sum(N_GWAS), + use_robust_var_est = use_robust_var_est) } Post_Model_Prob[C_id] <- temp_r2 Post_Model_Prob_Ratio[C_id] <- as.double(temp_r2 / Index_Model_Prob) @@ -416,7 +438,7 @@ mJAM_build_CS <- function(X_id, prev_X_list = NULL,All_id, PrCS_weights = "Pr(M_ mJAM_get_PrM_Wald(GItGI = GItGI_curr, GIty = GIty_curr, yty = yty_curr, yty_med = yty_med, rare_SNPs = rare_SNPs, N_GWAS = N_GWAS, C_id = X_id, prev_X_list = prev_X_list, - g = sum(N_GWAS)) + g = sum(N_GWAS), use_robust_var_est = use_robust_var_est) Std_CS_Prob[X_id] <- Post_Model_Prob[X_id]/sum_Post_CS_Porb } if(PrCS_weights == "Pr(M_C)"){ @@ -425,12 +447,13 @@ mJAM_build_CS <- function(X_id, prev_X_list = NULL,All_id, PrCS_weights = "Pr(M_ mJAM_get_PrM(GItGI = GItGI_curr, GIty = GIty_curr, yty = yty_curr, yty_med = yty_med, rare_SNPs = rare_SNPs, N_GWAS = N_GWAS, C_id = X_id, prev_X_list = prev_X_list, - g = sum(N_GWAS))$post_prob + g = sum(N_GWAS), use_robust_var_est = use_robust_var_est)$post_prob Std_CS_Prob[X_id] <- as.numeric(mJAM_get_PrM(GItGI = GItGI_curr, GIty = GIty_curr, yty = yty_curr, yty_med = yty_med,rare_SNPs = rare_SNPs, N_GWAS = N_GWAS, C_id = X_id, prev_X_list = prev_X_list, - g = sum(N_GWAS))$post_prob/sum_Post_CS_Porb) + g = sum(N_GWAS), + use_robust_var_est = use_robust_var_est)$post_prob/sum_Post_CS_Porb) } for(C_id in Testing_id){ @@ -439,13 +462,13 @@ mJAM_build_CS <- function(X_id, prev_X_list = NULL,All_id, PrCS_weights = "Pr(M_ temp_r2 <- mJAM_get_PrM_Wald(GItGI = GItGI_curr, GIty = GIty_curr, yty = yty_curr, yty_med = yty_med, rare_SNPs = rare_SNPs, N_GWAS = N_GWAS, C_id = C_id, prev_X_list = prev_X_list, - g = sum(N_GWAS)) + g = sum(N_GWAS), use_robust_var_est = use_robust_var_est) } if(PrCS_weights == "Pr(M_C)"){ temp_r2 <- mJAM_get_PrM(GItGI = GItGI_curr, GIty = GIty_curr, yty = yty_curr, yty_med = yty_med, rare_SNPs = rare_SNPs, N_GWAS = N_GWAS, C_id = C_id, prev_X_list = prev_X_list, - g = sum(N_GWAS))$post_prob + g = sum(N_GWAS), use_robust_var_est=use_robust_var_est)$post_prob } Post_Model_Prob[C_id] <- temp_r2 diff --git a/man/mJAM_build_CS.Rd b/man/mJAM_build_CS.Rd index 89958da..fbb2ff4 100644 --- a/man/mJAM_build_CS.Rd +++ b/man/mJAM_build_CS.Rd @@ -16,7 +16,8 @@ mJAM_build_CS( yty_med, N_GWAS, rare_SNPs = NULL, - Pr_Med_cut = 0.1 + Pr_Med_cut = 0.1, + use_robust_var_est = FALSE ) } \arguments{ @@ -43,6 +44,8 @@ mJAM_build_CS( \item{rare_SNPs}{A numeric vector of ID(s) for rare SNP(s) which we do not apply weighting. Instead, we use the individual estimate of yty for these SNPs for robustness.} \item{Pr_Med_cut}{The cutoff for Pr(Mediation); SNPs with Pr(Mediation) smaller than this cutoff will be assigned a Pr(CS) = 0 and thus not included in the credible set for the current index} + +\item{use_robust_var_est}{whether to use linear combination of median yty and individual yty.} } \value{ A table with the following columns: diff --git a/man/mJAM_get_PrM.Rd b/man/mJAM_get_PrM.Rd index 08fe21f..78bd91b 100644 --- a/man/mJAM_get_PrM.Rd +++ b/man/mJAM_get_PrM.Rd @@ -13,7 +13,8 @@ mJAM_get_PrM( C_id, prev_X_list = NULL, g = NULL, - rare_SNPs = NULL + rare_SNPs = NULL, + use_robust_var_est = FALSE ) } \arguments{ @@ -34,6 +35,8 @@ mJAM_get_PrM( \item{g}{The pre-specified `g` in `g`-prior formulation.} \item{rare_SNPs}{A numeric vector of ID(s) for rare SNP(s) which we do not apply weighting. Instead, we use the individual estimate of yty for these SNPs for robustness.} + +\item{use_robust_var_est}{whether to use linear combination of median yty and individual yty.} } \value{ \describe{ diff --git a/man/mJAM_get_PrM_Wald.Rd b/man/mJAM_get_PrM_Wald.Rd index b410541..e315a1e 100644 --- a/man/mJAM_get_PrM_Wald.Rd +++ b/man/mJAM_get_PrM_Wald.Rd @@ -13,7 +13,8 @@ mJAM_get_PrM_Wald( C_id, prev_X_list = NULL, g = NULL, - rare_SNPs = NULL + rare_SNPs = NULL, + use_robust_var_est = FALSE ) } \arguments{ @@ -34,6 +35,8 @@ mJAM_get_PrM_Wald( \item{g}{The pre-specified `g` in `g`-prior formulation.} \item{rare_SNPs}{A numeric vector of ID(s) for rare SNP(s) which we do not apply weighting. Instead, we use the individual estimate of yty for these SNPs for robustness.} + +\item{use_robust_var_est}{whether to use linear combination of median yty and individual yty.} } \value{ A numeric vector of posterior Pr(Model) for each SNPs in `C_id`.