R/uncor_CPBayes_functions.R

In CPBayes: Bayesian Meta Analysis for Studying Cross-Phenotype Genetic Associations

 ##=================================****** uncorrelated version of CPBayes ******================================##
 ## This function is the main MCMC function for implementing CPBayes in case of uncorrelated summary statistics.
 ## Uncorrelated summary statistics arise for separate case-control statistics without any overlapping subjects. 
 ## It calls the function: initiate_MCMC() from 'common_MCMC_functions.R' to initialize the parameters in the MCMC.
 ## It calls: uncorrelated_beta_update(), Z_update(), q_update(), de_update1(), de_update0() 
 ## from 'uncor_MCMC_functions.R' to update different parameters in MCMC.
 ## It calls select_subset(), overall_pleio_measure() from 'summary_functions.R' to summarize the MCMC data.
 ##=================================**********************************************===============================## 

 #library("MASS")

 CPBayes_uncor = function( variantName, traitNames, X, s.e., updateDE, MinSlabVar, MaxSlabVar, RP, burn.in )
 {
     ptm1 <- proc.time()
     set.seed(10)

     K = length(X)
     PPAj_thr = 0.20                                           ## PPAj threshold

     tau <- 0.01                                               ## choice of spike sd (var = tau^2)
     CentralSlabVar <- (MinSlabVar+MaxSlabVar)/2
     nonNullVar <- CentralSlabVar                              ## central choice of slab variance 
     de <- sqrt(tau^2/nonNullVar)                              ## 1/de = ratio of slab sd and spike sd

     ## v0 (minimum of slab variance - min_var), v1 (maximum of slab variance - max_var)
     min_var <- MinSlabVar
     max_var <- MaxSlabVar
     max_de <- tau/sqrt(min_var)                               ## maximum value of 'de'
     min_de <- tau/sqrt(max_var)                               ## mimimum value of 'de'
     shape1_de <- 1                                            ## shape1 parameter of the Beta prior of 'de' (shape2 parameter = 1, always)

    
     ## an informed initialization of the MCMC parameters 
     initiate <- initiate_MCMC( K, X, s.e. )
     beta <- initiate$beta
     Z <- initiate$Z
     q <- initiate$q
     nA <- initiate$K1_FDR                                     ## number of associated traits

     shape2 <- 1                                               ## shape2 parameter for Beta prior of q
     LB <- 0.1; UB <- 0.5;
     qm <- nA/K
     if(qm < LB) qm <- LB
     if(qm > UB) qm <- UB

     #qm <- 0.25
     shape1 <- (qm/(1-qm)) * shape2

     thinning <- 1                                             ## thinning period in the MCMC

     mcmc.samplesize <- (RP-burn.in) %/% thinning; Z.data <- matrix(0,mcmc.samplesize,K); row <- 0 ;
     sim.beta <- matrix(0,mcmc.samplesize,K); sample_probZ_zero <- matrix(0,mcmc.samplesize,K);

     for( rp in 1:RP )
     {
         ## Update beta
         beta = uncorrelated_beta_update(K, s.e., tau, de, Z, X)

         ## Update Z using the q-included version
         res_Z <- Z_update(K, q, tau, de, beta)
         
         ## using the q-integrated out version                                                                                                      
         #res_Z <- Z_integrated_update(K, log_ratio_p1, tau_const, de, beta)     

         ## collecting the otput from Z-update function
         Z <- res_Z$Z
         probZ_zero <- res_Z$prob

         q = q_update(K, Z, shape1, shape2)                    ## Update q
         ## Update de
         if(updateDE == TRUE){
             if(sum(Z) > 0) 
               de <- de_update1(min_de, max_de, shape1_de, beta, Z, tau) 
             else de <- de_update0(min_de, max_de, shape1_de)
         }

         ## collecting the MCMC sample obtained after the burn in period
         if(rp > burn.in && rp%%thinning == 0)
         {   
             row = row + 1 ;
             Z.data[row,] = Z ;
             sim.beta[row,] = beta ;
             sample_probZ_zero[row,] = probZ_zero;
         }
     }                                                        ## closing the loop for MCMC iterations     

     ##----------------------- Compute the summary obtained from the MCMC data ------------------------------------##
     ## selection of subset
     uncor.subset <- select_subset( K, Z.data, mcmc.samplesize )
     selected_traits <- NULL
     if( length(uncor.subset) > 0 ) 
       selected_traits <- traitNames[uncor.subset]
     
     ## extracting traits having PPAj > PPAj_thr         
     asso.pr = colSums(Z.data)/mcmc.samplesize
     which_traits = which(asso.pr > PPAj_thr)
     
     imp_PPAj = 0; imp_traits = 0; important_phenos = NULL
     
     if(length(which_traits) > 0){
     	imp_PPAj = asso.pr[which_traits]
     	imp_traits = traitNames[which_traits]
     	important_phenos = data.frame( traits = imp_traits, PPAj = imp_PPAj, stringsAsFactors = FALSE)
     } 
          
     ## calculate the Bayes factor and PPNA
     pleio_evidence <- overall_pleio_measure( K, shape1, shape2, sample_probZ_zero )
     log10_BF_uncor <- pleio_evidence$log10_BF
     PPNA.uncor <- pleio_evidence$PPNA    
     
     ptm2 <- proc.time()                                       ## time taken for the analysis
     ptm <- ptm2-ptm1
     #cat(" run time (in seconds):", "\n")
     #print(ptm2-ptm1)

     ## return the outputs. A post summary from the MCMC data can be computed for interesting variants
     data = list( variantName = variantName, log10_BF = log10_BF_uncor, locFDR = PPNA.uncor, 
            subset = selected_traits, important_traits = important_phenos, auxi_data = list( traitNames = traitNames, 
            K = K, mcmc.samplesize = mcmc.samplesize, PPAj = asso.pr, Z.data = Z.data, sim.beta = sim.beta, betahat = X, se = s.e.), runtime = ptm )

     return(data)
 }