# Dirac Spike Run a Gibbs sampler for a multivariate Bayesian sparse group selection model with Dirac spike prior for detecting pleiotropic effects on the traits. This function is designed for summary statistics containing estimated regression coefficients and their estimated covariance matrices. ```r DS( Betah, Sigmah, kappa0, sigma20, m, K, niter = 1000, burnin = 500, nthin = 2, nchains = 2, a1 = 0.1, a2 = 0.1, d1 = 0.1, d2 = 0.1, snpnames, genename ) ``` ## Arguments - `Betah`: A list containing m-dimensional vectors of the regression coefficients for K studies. - `Sigmah`: A list containing the positive definite covariance matrices (m*m-dimensional) which is the estimated covariance matrices of K studies. - `kappa0`: Initial value for kappa (its dimension is equal to nchains). - `sigma20`: Initial value for sigma2 (its dimension is equal to nchains). - `m`: Number of variables in the group. - `K`: Number of traits. - `niter`: Number of iterations for the Gibbs sampler. - `burnin`: Number of burn-in iterations. - `nthin`: The lag of the iterations used for the posterior analysis is defined (or thinning rate). - `nchains`: Number of Markov chains, when nchains>1, the function calculates the Gelman-Rubin convergence statistic, as modified by Brooks and Gelman (1998). - `a1, a2`: Hyperparameters of kappa. Default is a1=0.1 and a2=0.1. - `d1, d2`: Hyperparameters of sigma2. Default is d1=0.1 and d2=0.1. - `snpnames`: Names of variables for the group. - `genename`: Name of group. ## Returns * mcmcchain: The list of simulation output for all parameters. * Summary: Summary statistics for regression coefficients in each study. * Criteria: genename, snpnames, PPA, log10BF, lBFDR, theta. * Indicator: A table containing m rows of binary indicators for each study, the number of studies with nonzero signal and having pleiotropic effect by credible interval (CI). The first K columns show nonzero signals, K+1 th column includes the number of studies with nonzero signal and the last column shows an indicator for having pleiotropic effect of each SNP. ## Details Let betah_k, k=1,...,K be a m-dimensional vector of the regression coefficients for the kth study and Sigmah_k be its estimated covariance matrix. The hierarchical set-up of DS prior, by considering summary statistics (betah_k and Sigmah_k, k=1,...,K) as the input of the method, is given by: betah _k ~ (1 - kappa) delta_0(betah_k) + kappa N_m(0,sigma2 I_m ), kappa ~ Beta(a_1,a_2), sigma2 ~ inverseGamma (d_1,d_2). where delta_0(betah_k) denotes a point mass at 0, such that delta_0(betah_k)=1 if beta_k=0 and delta_0(betah_k)=0 if at least one of the $m$ components of beta_k is non-zero. ## Examples ```r ############################# Gene DNAJC1 ############################################### data(DNAJC1) Breast <- DNAJC1$Breast Thyroid <- DNAJC1$Thyroid genename <- "DNAJC1" snpnames <- Breast$snp Betah <- list(Breast$beta, Thyroid$beta) Sigmah <- list(diag(Breast$se^2), diag(Thyroid$se^2)) K <- 2 m <- 14 RES <- DS(Betah, Sigmah, kappa0 = 0.5, sigma20 = 1, m = m, K = K, niter = 2000, burnin = 1000, nthin = 2, nchains = 1, a1 = 0.1, a2 = 0.1, d1 = 0.1, d2 = 0.1, snpnames, genename ) ## Not run: print(RES) RES1 <- DS(Betah, Sigmah, kappa0 = c(0.2, 0.5), sigma20 = c(1, 2), m = m, K = K, niter = 2000, burnin = 1000, nthin = 2, nchains = 2, a1 = 0.1, a2 = 0.1, d1 = 0.1, d2 = 0.1, snpnames, genename ) print(RES1) ################### Simulated summary level data with K=5 ############################### data(Simulated_summary) genename <- Simulated_summary$genename snpnames <- Simulated_summary$snpnames Betah <- Simulated_summary$simBeta Sigmah <- Simulated_summary$simSIGMA K <- 5 m <- 10 RES <- DS(Betah, Sigmah, kappa0 = 0.5, sigma20 = 1, m = m, K = K, niter = 2000, burnin = 1000, nthin = 2, nchains = 1, a1 = 0.1, a2 = 0.1, d1 = 0.1, d2 = 0.1, snpnames, genename ) print(RES) RES1 <- DS(Betah, Sigmah, kappa0 = c(0.2, 0.5), sigma20 = c(1, 2), m = m, K = K, niter = 2000, burnin = 1000, nthin = 2, nchains = 2, a1 = 0.1, a2 = 0.1, d1 = 0.1, d2 = 0.1, snpnames, genename ) print(RES1) ################################### Gene PARP2 ########################################## library(BhGLM) data(PARP2) Breast <- PARP2$Breast Thyroid <- PARP2$Thyroid genename <- "PARP2" snpnames <- c("rs3093872", "rs3093921", "rs1713411", "rs3093926", "rs3093930", "rs878156") Fit1 <- BhGLM::bglm(y1 ~ ., family = binomial(link = "logit"), data = Breast) Betah1 <- Fit1$coefficients[-1] Sigmah1 <- cov(coef(arm::sim(Fit1)))[-1, -1] Fit2 <- BhGLM::bglm(y2 ~ ., family = binomial(link = "logit"), data = Thyroid) Betah2 <- Fit2$coefficients[-1] Sigmah2 <- cov(coef(arm::sim(Fit2)))[-1, -1] Betah <- list(Betah1, Betah2) Sigmah <- list(Sigmah1, Sigmah2) K <- 2 m <- 6 RES <- DS(Betah, Sigmah, kappa0 = 0.5, sigma20 = 1, m = m, K = K, niter = 2000, burnin = 1000, nthin = 2, nchains = 1, a1 = 0.1, a2 = 0.1, d1 = 0.1, d2 = 0.1, snpnames, genename ) print(RES) RES1 <- DS(Betah, Sigmah, kappa0 = c(0.2, 0.5), sigma20 = c(1, 2), m = m, K = K, niter = 2000, burnin = 1000, nthin = 2, nchains = 2, a1 = 0.1, a2 = 0.1, d1 = 0.1, d2 = 0.1, snpnames, genename ) print(RES1) ########### Simulated individual level data with K=3 and continuous phynotype ########### library(BhGLM) data(Simulated_individual) Study1 <- Simulated_individual$Study1 Study2 <- Simulated_individual$Study2 Study3 <- Simulated_individual$Study3 K <- 3 m <- 30 genename <- "Simulated" snpnames <- sprintf("SNP%s", seq(1:m)) Fit1 <- BhGLM::bglm(Y1 ~ ., family = gaussian, data = data.frame(Study1)) Betah1 <- Fit1$coefficients[-1] Sigmah1 <- cov(coef(arm::sim(Fit1)))[-1, -1] Fit2 <- BhGLM::bglm(Y2 ~ ., family = gaussian, data = data.frame(Study2)) Betah2 <- Fit2$coefficients[-1] Sigmah2 <- cov(coef(arm::sim(Fit2)))[-1, -1] Fit3 <- BhGLM::bglm(Y3 ~ ., family = gaussian, data = data.frame(Study3)) Betah3 <- Fit3$coefficients[-1] Sigmah3 <- cov(coef(arm::sim(Fit3)))[-1, -1] Betah <- list(Betah1, Betah2, Betah3) Sigmah <- list(Sigmah1, Sigmah2, Sigmah3) RES <- DS(Betah, Sigmah, kappa0 = 0.5, sigma20 = 1, m = m, K = K, niter = 2000, burnin = 1000, nthin = 2, nchains = 1, a1 = 0.1, a2 = 0.1, d1 = 0.1, d2 = 0.1, snpnames, genename ) print(RES) RES1 <- DS(Betah, Sigmah, kappa0 = c(0.2, 0.5), sigma20 = c(1, 2), m = m, K = K, niter = 2000, burnin = 1000, nthin = 2, nchains = 2, a1 = 0.1, a2 = 0.1, d1 = 0.1, d2 = 0.1, snpnames, genename ) print(RES1) ########### Simulated individual level data with K=2 and gene expression data ########### library(BhGLM) data(Simulated_individual_survival) Study1 <- Simulated_individual_survival$Study1 Study2 <- Simulated_individual_survival$Study2 K <- 2 m <- 10 genename <- "Simulated" snpnames <- sprintf("G%s", seq(1:m)) Fit1 <- BhGLM::bcoxph(Study1$T ~ Study1$X) Betah1 <- Fit1$coefficients Sigmah1 <- Fit1$var Fit2 <- BhGLM::bcoxph(Study2$T ~ Study2$X) Betah2 <- Fit2$coefficients Sigmah2 <- Fit2$var Betah <- list(Betah1, Betah2) Sigmah <- list(Sigmah1, Sigmah2) RES1 <- DS(Betah, Sigmah, kappa0 = c(0.2, 0.5), sigma20 = c(1, 2), m = m, K = K, niter = 2000, burnin = 1000, nthin = 2, nchains = 2, a1 = 0.1, a2 = 0.1, d1 = 0.1, d2 = 0.1, snpnames, genename ) print(RES1) ## End(Not run) ``` ## References 1. Baghfalaki, T., Sugier, P. E., Truong, T., Pettitt, A. N., Mengersen, K., & Liquet, B. (2021). Bayesian meta analysis models for cross cancer genomic investigation of pleiotropic effects using group structure. **Statistics in Medicine**, 40 (6), 1498-1518. ## Author(s) Taban Baghfalaki.

DS function