Skip to content

Commit

Permalink
added use_robust_var_est option in mJAM_build_CS
Browse files Browse the repository at this point in the history
  • Loading branch information
@jiayi-s
jiayi-s committed Nov 30, 2023
1 parent 4cec87e commit 47b00c9
Show file tree
Hide file tree
Showing 5 changed files with 53 additions and 20 deletions.
3 changes: 2 additions & 1 deletion R/mJAM_Forward.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
55 changes: 39 additions & 16 deletions R/mJAM_build_CS.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)}
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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

Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -360,15 +378,17 @@ 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
}
if(PrCS_weights == "Pr(Wald)"){
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
}

Expand All @@ -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)
Expand All @@ -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)"){
Expand All @@ -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){
Expand All @@ -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
Expand Down
5 changes: 4 additions & 1 deletion man/mJAM_build_CS.Rd
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

5 changes: 4 additions & 1 deletion man/mJAM_get_PrM.Rd
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

5 changes: 4 additions & 1 deletion man/mJAM_get_PrM_Wald.Rd
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

0 comments on commit 47b00c9

Please sign in to comment.